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Diversity  is  an  average  property  of  a  community 
and  is  identified  as  the  species  variety  and 
abundance.  The  study  of  biodiversity  is 
important  because  it  is  one  of  the  central 
themes  in  ecology;  the  diversity  of  a  system  is 
often  seen  as  an  indicator  of  the  well-being  of 
the  system.  In  this  study,  we  selected  and 
defined  theoretical  statistical  measures  of  plant 
diversity  and  developed  theoretical  dynamics 
models  for  plant  communities.  This  project 
provides  a  new  approach  to  measuring  plant 
diversity  through  time. 

The  stochastic  dynamics  modeled  in  this  project 
contain  two  main  components:  deterministic 
processes  and  stochastic  processes. 


Stochastic  dynamics  models  make  it  possible  to 
simulate  plant  diversity  changes  through  time 
even  without  long-term  observed  data.  Because 
both  the  time  and  space  dimensions  are 
included  in  these  stochastic  dynamics  models, 
they  have  more  extensive  uses  in  assessing 
and  monitoring  plant  communities.  The 
potential  applications  of  these  models  include 
(1)  providing  standard  diversity  measures,  (2) 
monitoring  the  development  of  plant 
communities  in  terms  of  species  diversity  and 
structure  diversity,  (3)  testing  the  significance  of 
the  influence  of  human  activities  on  plant 
communities,  and  (4)  estimating  the 
rehabilitation  time  of  a  disturbed  plant 
community. 


Approved  for  public  releeise:  distribution  is  unlimited.  W(WW.cecer.army.mil/T echReports 

W3C  QOAiiii  X  3 


CERL  TR  00-5 


Foreword 


This  study  was  conducted  for  the  Office  of  the  Directorate  of  Environmental  Pro¬ 
gram  (DAIM),  Assistant  Chief  of  Staff  (Installation  Management)  (ACS(IM))  im- 
der  Project  4A1622720A896,  “Environmental  Quality  Technology;”  Work  Unit 
CN-T09,  “Installation  Capacity  Factors.”  The  technical  monitor  was  Dr.  Victor 
E.  Diersing,  DAIM-ED-N. 

The  work  was  performed  by  the  Ecological  Processes  Branch  (CN-N)  of  the  In¬ 
stallations  Division  (CN),  U.S.  Army  Construction  Engineering  Research  Labo¬ 
ratory  (CERL).  The  CERL  Principal  Investigator  was  Alan  B.  Anderson.  Steve 
Hodapp  is  Chief,  CECER-CN-N,  and  Dr.  John  Bandy  is  Chief,  CECER-CN.  The 
technical  editor  was  Gloria  Wienke,  Information  Technology  Laboratory. 

Part  of  this  work  was  done  by  Xianchi  Cao  and  Dr.  George  Gertner  of  the  Uni¬ 
versity  of  Illinois  Department  of  Natural  Resources  and  Environmental  Sciences. 
Bruce  MacAllister  worked  on  the  project  as  an  Oakridge  Associated  Universities 
post-graduate  research  fellow.  He  also  provided  the  cover  graphic  for  the  report. 

The  authors  would  like  to  thank  White  Sands  Missile  Range  ITAM  personnel  for 
access  to  LCTA  data  and  review  of  the  report. 

The  Acting  Director  of  CERL  is  Dr.  Alan  W.  Moore. 

CERL  is  an  element  of  the  U.S.  Army  Engineer  Research  and  Development  Cen¬ 
ter  (ERDC),  U.S.  Army  Corps  of  Engineers.  The  Acting  Director  of  ERDC  is  Dr. 
Lewis  E.  Link  and  the  Commander  is  COL  Robin  R.  Cababa,  EN. 


DISCLAIMER 

The  contents  of  this  report  are  not  to  be  used  for  advertising,  publication,  or  promotional  purposes.  Citation  of  trade  names 
does  not  constitute  an  official  endorsement  or  approval  of  the  use  of  such  commercial  products.  All  product  names  and 
trademarks  cited  are  the  property  of  their  respective  ovwiers. 

The  findings  of  this  report  are  not  to  be  construed  as  an  official  Department  of  the  Army  position  unless  so  designated  by 
other  authorized  documents. 


DESTROY  THIS  REPORT  WHEN  IT  IS  NO  LONGER  NEEDED.  DO  NOT  RETURN  IT  TO  THE  ORIGINATOR. 


CERL  TR  00-5 


3 


Contents 

Foreword . 2 

List  of  Figures  and  Tables . 1 . 8 

1  Introduction . 9 

Background . 9 

Objectives . 9 

Approach . "lO 

Scope . 19 

Mode  of  Technology  Transfer . 10 

2  Ecological  Diversity  and  Its  Measurement . 11 

The  Concept  of  Diversity  and  its  Measurements . 11 

Indices  of  Species  Richness . 12 

Species  Abundance  Models . 12 

Diversity  Indices. . 1^ 

New  Models  of  Diversity . 16 

3  Stochastic  Dynamic  Models . 19 

A  Simulation  Model  of  Population  Dynamics . ••••20 

A  Discrete  Model  of  Plant  Population  Dynamics . 20 

The  Continuous  Model  of  Plant  Population  Dynamics . 22 

Demographic  Stochasticity . 23 

The  Poisson  Model  of  Immigration . 23 

The  Simple  Birth  and  Death  Process . 24 

Environmental  Stochasticity . 27 

Stochastic  Differential  Equations  (1) . 25 

Stochastic  Differential  Equations  (2) . 34 

Joint  Effect  of  Demographic  and  Environmental  Stochasticity . 37 

4  Stochastic  Estimation  and  Identification . 40 

Optimal  Estimation  . . 41 

Optimal  Linear  Filtering . 42 

Optimal  Prediction . 46 

Nonlinear  Estimation . 47 

Stochastic  Identification . 50 


Least  Square  Mean  Estimates  (LSE) . 50 

Conditional  Mean  Estimates . . 52 

5  STELLA  II  Modeling  Process . 55 

The  STELLA  Model.. . 57 

Setting  of  the  STELLA  model . 60 

Giobal  Variables . 60 

Local  Variables . 61 

Kalman  filter. . 61 

Outputs . 62 

6  A  Case  Study:  White  Sands  Instaiiation . 63 

Classification  of  Plant  Communities . 63 

Standardization  of  Data  by  Quadrat . 64 

Clustering  Method . 65 

Average  Linkage . 66 

Centroid  Method . 66 

Ward’s  Minimum-variance  Method . 67 

Classification  Results . 67 

Simulation  Results . 71 

7  Summary . 80 

References . 82 

DISTRIBUTION . 86 

Report  Documentation  Page . 87 


CERL  TR  00-5  _ 5 


List  of  Figures  and  Tables 

Figures 

1  Sequence  of  choosing  a  stochastic  model . 55 

2  The  STELLA  map  of  flow  structure . 58 

3  The  STELLA  structure  of  the  Kalman  filter. . 59 

4  Simulation  results  for  plant  community  type  1 . 71 

5  Simulation  results  for  plant  community  type  2 . 71 

6  Simulation  results  for  plant  community  type  3 . 72 

7  Simulation  results  for  plant  community  type  4 . 72 

8  Simulation  results  for  plant  community  type  5 . 73 

9  Simulation  results  for  plant  community  type  6 . 73 

1 0  Simulation  results  for  plant  community  type  7 . 74 

11  Simulation  results  for  plant  community  type  8 . = . 74 

1 2  Simulation  results  for  plant  community  type  9 . 75 

1 3  Simulation  results  for  plant  community  type  1 0 . 75 

1 4  Simulation  results  for  plant  community  type  11 . 76 

1 5  Simulation  results  for  plant  community  type  12 . 76 

1 6  Simulation  results  for  plant  community  type  1 3 . 77 

1 7  Simulation  results  for  plant  community  type  14 . 77 

1 8  Simulation  results  for  plant  community  type  15 . 78 

1 9  Simulation  results  for  plant  community  type  16 . ; . 78 

20  Simulation  results  for  plant  community  type  1 7 . 79 

Tables 

1  Comparison  of  cluster  analysis  results  for  1 989 . 68 

2  Comparison  of  cluster  analysis  results  for  1 992 . 69 

3  Vegetation  types  identified  at  White  Sands  Missile  Range . 70 


CERL  TR  00-5 


7 


1  Introduction 

Background 

The  United  States  Army  is  responsible  for  managing  over  12  million  acres  of 
land.  The  Army’s  land  management  objective  is  to  maintain  realistic  military 
training  and  testing  environments  while  promoting  land  stewardship.  To  ac- 
comphsh  this  objective,  the  U.S.  Army  Land  Condition  Trend  Analysis  (LCTA) 
program  was  developed  at  the  U.S.  Army  Construction  Engineering  Research 
Laboratory  (CERL)  under  the  sponsorship  of  the  former  U.S.  Army  Engineering 
and  Housing  Support  Center  (USAEHSC)  as  a  means  to  inventory  and  monitor 
natural  resources  on  mihtary  installations.  LCTA  uses  standard  methods  to 
collect,  analyze,  and  report  natural  resources  data  (Diersing,  Shaw,  and  Tazik 
1992)  and  is  the  Army’s  standard  for  land  inventory  and  monitoring  (Technical 
Note  [TN]  420-74-3).  Over  50  military  installations  and  training  areas  in  the 
United  States  and  Germany  have  begun  or  plan  to  implement  LCTA.  LCTA  data 
sets  currently  exist  for  more  than  40  installations  and  contain  1  to  10  years  of 
monitoring  data.  Lands  inventoried  as  part  of  the  LCTA  program  include  Army 
Materiel  Command  (AMC),  Forces  Command  (FORSCOM),  Training  and  Doc¬ 
trine  Command  (TRADOC),  and  National  Guard  Bureau  installations.  Over  75 
percent  of  the  Army’s  land  base  is  represented  by  LCTA  data  (Shaw  and  Kowal¬ 
ski  1996). 

An  informal  review  of  installation  ITAM  personnel  indicated  an  interest  in  esti¬ 
mating  plant  diversity  using  LCTA  data  and  modehng  changes  in  plant  diversity 
that  result  from  alternative  land  uses. 


Objectives 

The  objective  of  this  project  was  to  develop  and  test  methodology  to  model 
changes  in  plant  diversity  using  standard  data  from  the  U.S.  Armys  Land  Con¬ 
dition  Trend  Analysis  (LCTA)  program.  Specifically,  stochastic  models  (those  in¬ 
volving  random  variables)  of  plant  diversity  were  to  be  developed  using  data 
from  the  White  Sands  Missile  Range,  New  Mexico,  LCTA  program. 
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Approach 

A  literature  survey  was  conducted  to  identify  methods  of  characterizing  and 
modeling  plant  diversity.  Based  on  results  of  the  literatvire  review,  modehng 
techniques  were  selected  to  model  plant  diversity.  Land  Condition  Trend  Analy¬ 
sis  data  from  the  White  Sands  Missile  Range,  New  Mexico,  was  then  used  to  de¬ 
velop  stochastic  plant  diversity  models  for  selected  plant  commxinities. 


Scope 

The  modeling  techniques  described  in  this  report  are  applicable  to  any  military 
installation  collecting  LCTA  data.  The  specific  model  documented  is  only  appli¬ 
cable  to  White  Sands  Missile  Range. 


Mode  of  Technology  Transfer 

The  information  in  this  report  will  be  provided  to  Army  ITAM  personnel  respon¬ 
sible  for  using  LCTA  data.  The  information  will  also  be  provided  to  organiza¬ 
tions  responsible  for  developing  and  refining  natural  resoxirces  conservation 
methodologies  through  hard  copy  reports  and  through  the  CERL  web  site. 

This  study  is  part  of  a  larger  research  effort  to  develop  and  field  LCTA-related 
data  apphcations.  Models  and  analysis  techniques  described  in  this  report  are 
being  incorporated  into  other  modeling  efforts  and  systems. 


CERL  TR  00-5 


9 


2  Ecological  Diversity  and  its 
Measurement 

Although  the  study  of  ecological  diversity  is  only  about  30  years  old,  it  has  been 
extensively  represented  in  the  literature  during  that  time.  There  are  three  rea¬ 
sons  ecologists  are  so  interested  in  this  topic.  First,  diversity  is  one  of  the  cen¬ 
tral  themes  in  ecology.  Ecology  is  the  scientific  study  of  the  distribution  and 
abundance  of  organisms  and  the  interrelationship  between  the  orgeinisms  and 
their  surroimdings.  Diversity  studies,  which  address  the  variety  and  abimdance 
of  organisms,  represent  a  major  field  of  study  in  ecology.  Second,  diversity  is  of¬ 
ten  seen  as  an  indicator  of  the  health  of  an  ecosystem.  Studies  have  shown  that 
pollution  and  disturbance  reduce  the  richness  and  variety  of  the  natural  ecologi¬ 
cal  communities.  The  loss  of  natural  habitat  and  species  extinction  aroimd  the 
world  have  served  to  focus  international  attention  on  the  issue  of  diversity. 
Third,  considerable  debate  surroxmds  the  measurement  of  diversity.  On  the  sur¬ 
face,  biodiversity  seems  to  be  a  straightforward  concept.  Most  people  have  an 
intmtive  sense  of  the  word.  They  wotild  acknowledge  that  tropical  rain  forests 
harbor  more  species  than  temperate  woodlands  and  are  therefore  more  biologi¬ 
cally  diverse.  However,  the  more  we  look  at  diversity,  the  less  clearly  defined  it 
appears  to  be  because  diversity  can  be  measxired  in  so  many  different  ways.  A 
more  in-depth  study  of  diversity  could  reveal  new  and  unexpected  relationships 
between  species  and  ultimately  lead  to  a  better  understanding  of  the  mecha¬ 
nisms  involved.  The  study  of  ecological  diversity  over  the  past  30  years  has 
raised  three  main  questions:  what  is  diversity,  how  is  it  measiired,  and  what  is 
its  value  in  practice? 


The  Concept  of  Diversity  and  its  Measurements 

Diversity  is  one  property  of  a  biological  community  and  consists  of  two  compo¬ 
nents:  variety  and  abimdance.  A  large  number  of  diversity  measures  have  been 
devised  by  interpreting  the  relationship  between  variety  and  abundance  in  dif¬ 
ferent  ways.  Magurran  (1988)  divides  the  measurement  of  species  diversity  into 
three  main  categories.  First  are  the  species  richness  indices.  These  indices  are 
essentially  measures  of  the  number  of  species  in  an  ecosystem.  The  second  cate¬ 
gory  of  diversity  measures  includes  those  models  that  describe  the  distribution  of 
species  abimdance.  The  third  category  is  the  diversity  indices  based  on  the  rela- 
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tive  species  abundance.  These  indices,  like  the  Shannon  index  and  the  Simpson 
index,  consolidate  species  richness  and  evenness  into  a  single  figure. 


Indices  of  Species  Richness 


Species  richness  is  defined  as  the  number  of  species  or  species  density  in  the 
community.  If  a  complete  catalogue  of  species  in  a  community  can  be  obtained  (it 
is  possible  for  small  communities),  the  number  of  species  provides  some  measme 
of  understanding  diversity.  Because  most  natural  communities  are  very  large, 
however,  only  a  limited  number  of  species  can  be  counted  and  identified  by  sam¬ 
pling.  Therefore,  species  density  (defined  as  the  munber  of  species  per  unit  area, 
biomass,  or  number  of  individuals),  is  commonly  used  as  species  richness  (Hml- 
bert  1971;  Homer  1976;  Kempton  and  Wedderbum  1979;  Kershaw  and  Loony 
1985).  Species  density  estimated  by  samphng  varies  with  sample  size  and  sam¬ 
ple  distribution.  To  cope  with  this  problem  Sander  (1968)  devised  the  techmque 
of  Rarefaction  for  calculating  the  number  of  species  expected  for  all  samples  with 
standard  sample  size.  Hurlbert  (1971)  improved  Sanders’  Rarefaction  to  produce 
an  unbiased  estimate  of  the  number  of  species.  Instead  of  using  the  number  of 
species  or  the  species  density,  others  have  used  some  simple  indices  derived  from 
a  combination  of  the  number  of  species  (S)  and  the  number  of  total  individuals 
(N)  to  represent  species  richness.  Such  indices  include  Margalefis  index  (Clifford 
and  Stephenson  1975)  and  Menhinck’s  index  (Whittaker  1970).  The  Margalef 
index  (Dj,g)  and  Menhinck  index  (D^,^)  are,  respectively,  defined  by: 


(la) 


D 


Mg  ” 


s-\ 

InN 


(lb) 


D. 


S 

^fN 


Species  Abundance  Models 

Species  richness  may  be  intmtive  and  easy  to  calculate  but  it  does  not  contain 
any  information  of  the  relative  abundance  or  distribution  of  species.  In  fact,  spe¬ 
cies  distributions  are  often  more  meaningful  in  explaining  natural  communities. 
Kempton  and  Wedderbum  (1979)  observe  that  a  distribution  of  species  is  often  a 
more  sensitive  measure  of  environmental  disturbance  than  species  richness 
alone.  In  an  early  study,  Fisher,  Corbet,  and  Williams  (1943)  fotmd  that  patterns 
occur  in  species  distribution.  It  is  very  rare  to  have  an  equal  number  of  indi¬ 
viduals  for  £ill  species.  Instead,  a  few  species  would  be  very  abimdant,  some 
wotdd  have  a  medium  abundance,  while  most  species  have  only  a  few  individu¬ 
als.  This  observation  led  to  the  development  of  species  abvmdance  models. 
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Pielou  (1975)  developed  six  models  of  species  abiindance.  They  are  the  niche 
preemption  model,  broken  stick  model,  overlapping  niche  model,  truncated  nega¬ 
tive  binomial  distribution,  log-series  distribution,  and  log-normal  distribution. 
In  a  monograph  of  stochastic  abundance  models,  Engen  (1978)  lists  nine  mathe¬ 
matical  distributions  related  to  species  abundance.  These  are  the  gamma  distri¬ 
bution,  first  kind  beta  distribution,  second  kind  beta  distribution,  log-normal  dis¬ 
tribution,  Dirichlet  distribution,  negative  binomial  distribution,  logarithmic 
series  distribution,  negative  binomial  beta  distribution,  and  Poisson  log-normal 
distribution.  A  niimber  of  other  species  abimdance  models  also  appear  in  litera¬ 
ture  (Zipf  1965;  Mandelbrot  1977;  Gray  1988).  In  practice,  four  main  models 
(Magurran  1988)  characterize  species  diversity:  the  log-normal  distribution, 
geometric  series,  logarithmic  series,  and  MacArthur’s  broken  stick  model. 

The  geometric  series  is  based  on  the  hypothesis  that  every  species  ranked  from 
the  most  to  the  least  abimdant  take  the  same  proportion  (k)  of  the  remainder. 
The  ranked  abvmdance  list  is  k,  k(l-k),  k(l-k)^  ...  Ml-k)®"^,  (l-k)®^\  May  (1975) 
gave  the  probabihty  distribution  of  such  ranked  list,  F(x),  which  is  defined  as  the 
probability  that  a  randomly  chosen  species  has  size  less  than  x. 


(2) 


F(x)  =  C- 


ln(x) 

ln(l-A:) 


C  is  a  constant.  Geometric  series  pattern  of  species  abtmdance  is  fotmd  primar¬ 
ily  in  species-poor  environments  or  in  the  very  early  stages  of  a  succession  (Whit¬ 
taker  1965, 1970,  1972). 

Fisher,  Corbet,  and  Williams  (1943)  derived  a  log-series  model  to  describe  the 
species  abundance  of  Malayan  Lepidoptera.  This  log-series  model  represented 
the  first  attempt  to  describe  mathematically  the  relationship  between  the  num¬ 
ber  of  species  and  the  number  of  individuals  in  those  species.  In  the  log-series, 
the  expected  frequency  of  a  species  with  abxmdance  x  is  given  by: 

(3)  fx  =  ,  for  X  =  1,  2, ... 

X 

The  variable  b  (0  <  b  <  1)  is  a  constant  that  is  dependent  on  the  sample  size,  and 
a  (a  >  0)  is  a  constant  determined  by  the  characteristics  of  the  community.  When 
the  abundance  of  each  species  is  plotted  on  a  logarithmic  scale  in  rank,  the  log- 
series  approximates  a  straight  line  with  a  slope  of  -a  (Taylor,  Kempton,  and 
Woiwod  1976).  The  log-series  provides  a  statistical  satisfactory  description  of 
samples  from  a  wide  range  of  communities  (WilHams  1964;  Kempton  and  Taylor 
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1974;  Gray  1978).  It  can  be  used  for  small,  stressed,  or  pioneer  communities 
(May  1975;  Whittaker  1965),  and  also  for  “neutral”  imstressed  conummities 
(Caswell  1976). 

Preston  (1948,  1962)  observed  that  in  many  samples  middle-ranked  species  were 
relatively  numerous,  and  there  were  fewer  rare  species  than  the  log-series  dis¬ 
tribution  predicted.  By  the  log  transformation  of  number  of  individuals,  Preston 
foimd  that  the  number  of  species  always  distributed  with  a  truncated  normal 
distribution.  In  a  log-normal  distribution,  the  probability  density  function  is: 


(4) 


f(x)  = 


x^lm^ 


exp 


-(Inx-m)" 

2^^ 


The  expected  size  of  species  is  calculated  by: 


(5) 


E(S)  =  exp 


<T 

m-{ — 
2 


The  log-normal  distribution  rises  from  the  statistical  properties  of  large  munbers 
and  as  a  consequence  of  the  Central  Limit  Theorem  (May  1975).  Thus,  it  is  con¬ 
sidered  to  be  the  descriptor  of  large  and  mature  natm*al  conummities  (May  1975; 
WThittaker  1972;  Gray  1978;  Preston  1980;  Magurran  1988). 


MacArthur  (1957)  first  proposed  the  broken  stick  model.  In  this  model,  the  re- 
soimce  is  likened  to  a  stick  broken  randomly  and  simultaneously  into  S  distinct 
segments.  The  lengths  of  the  segments  represent  the  “sizes”  of  the  species.  Ac¬ 
cording  to  the  model,  the  expected  size  of  the  i-th  species,  xi,  is: 

<6)  E(x,)=ii  I 

The  broken  stick  model  reflects  a  much  more  equitable  state  of  affairs  than  those 
suggested  by  the  log-normal,  log-series,  and  geometric  series.  It  has  good  fits  for 
the  conummities  with  a  few  species  and  relative  high  evenness  between  species 
(May  1974;  Pielou  1975). 

Diversity  indices 

The  third  kind  of  diversity  meastu’es  include  those  indices  based  on  the  propor¬ 
tional  abimdance  of  species.  Although  species  abvmdance  distributions  provide 
the  fullest  description  of  diversity  data,  there  are  times  when  a  single  diversity 
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index  is  required.  When  communities  do  not  fit  one  model  and  it  is  desired  to 
compare  them  by  means  of  diversity,  indices  based  on  the  proportional  ahrm- 
dance  of  species  provide  a  solution  to  this  problem.  The  most  commonly  used 
indices  are  the  Shannon  index  and  Simpson  index. 

The  Shannon  index  is  defined  as: 

(7)  h’=-2;  (p,tep,) 

The  parameter  p.  is  the  proportion  of  community  members  in  the  i-th  class. 
Shannon  originally  proposed  this  diversity  index  as  a  measure  of  the  information 
content  of  a  code.  The  Shannon  index  assumes  that  individuals  are  randomly 
sampled  from  an  “infinitely  large”  population  and  all  species  are  represented  in 
the  sample. 

Using  the  fact  that  the  probability  of  drawing  two  successive  individuals  be¬ 
longing  to  the  same  species  in  a  random  sampling  is  p**.,  Simpson  (1949)  sug¬ 
gested  a  statistic,  D,  that  has  the  form  of: 

(8) 

i=l 

The  parameter  Pj  is  the  proportion  of  individuals  in  the  i-th  species  and  S  is  the 
total  number  of  species  in  the  commvmity.  This  statistic  measures  a  property 
that  is  opposite  to  the  diversity.  The  diversity  index  corresponding  to  the  statis¬ 
tic  D,  the  Simpson  index,  is  then  given  by: 

(9)  H  =  -InD 

The  Shannon  index  and  Simpson  index  are  two  special  cases  of  a  more  general 
class  of  functions  (H)  used  in  mathematical  theory  of  information  (Pielou  1975). 

Inctpf) 

(10)  H,  =— a - 

l-a 

(11)  H,  =|ijjj(H„)  =  -^p,.ln(p,.)  (Shannon  index) 

a~>\  /=:1 

f  =  l 


(12) 


(Simpson  index) 
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Using  the  concept  of  rarity,  Dennis  and  Patil  (1979)  find  a  class  of  diversity  func¬ 
tions  that  also  leads  to  the  Shannon  index  and  Simpson  index. 


(13) 

Rarity  =  R(pi)  =  p"l[l-pf] 

(14) 

s 

Diversity  =  Ap  =  2]piR(Pi ) 

(15) 

A-i  =  Ap=.i  =  S-1 

(Species  richness  index) 

(16) 

^0  =  = 

2=1 

(Shannon  index) 

(17) 

A,  =A^.,  =1-XP/ 

(Simpson  index) 

/=! 


There  are  other  diversity  indices  such  as  the  McIntosh  index  and  Berger-Parker 
index.  McIntosh  (1967)  found  that  the  Euclidean  distance  of  the  assemblage 
from  the  origin  could  be  used  as  a  measure  of  diversity.  Berger  and  Parker 
(1970)  use  the  proportional  importance  of  the  most  abimdant  spices  as  the  diver¬ 
sity  measure. 


New  Models  of  Diversity 

Hughes  (1984,  1986)  has  completed  research  on  the  diversity  measures  with  a 
dynamics  model.  In  this  dynamics  model,  the  abundance  (n)  of  the  i-th  species  at 
time  t+1  is  calculated  from  the  expression: 

(18)  <1  =  S[ni  +  R(1  -I-  Zni  )^~^*] 

jK 

The  variable  S  is  the  survival  petrameter,  R  is  the  recruit  parameter,  Z  is  the  ag¬ 
gregation  parameter,  and  is  the  number  of  total  individuals.  This  d3mamics 
model  simulates  the  development  and  the  progression  of  a  theoretical  commu¬ 
nity  through  time.  It  defines  a  “community”  with  a  variable  number  of  species 
and  species  abundance. 

Diversity  itself  is  also  divided  into  two  types:  species  diversity  and  spatial  diver¬ 
sity  (structure  and  habitat  diversity).  The  diversity  discussed  above  is  mainly 
concerned  with  species  diversity.  Although  species  diversity  is  the  more  impor¬ 
tant,  structure  and  habitat  diversity  has  special  use  in  ecology.  Habitat  diversity 
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has  been  used  as  an  important  component  of  wildlife  conservation  evaluation 
(Fuller  and  Langslow  1986;  Pearsall,  Durham,  and  Eagar  1986;  Usher  1986). 
MacArthur  and  MacArthur  (1961)  found  that  the  structm-al  diversity  of  a  tem¬ 
perate  woodland  in  North  America  was  a  much  better  predictor  of  bird  diversity 
than  the  plant  species  diversity.  Southwood,  Brown,  and  Reader  (1979)  reported 
finding  a  close  relationship  between  plant  spatial  diversity  and  insect  diversity 
in  woodland  succession.  Elton  and  Liller  (1954)  divided  habitat  diversity  into 
four  levels.  The  first  level  is  the  major  habitat  system  (e.g.,  terrestrial  or 
aquatic).  The  major  habitat  system  is  then  divided  into  formation  type  (e.g., 
woodland  or  open  land).  When  categorizing  the  formation  type,  the  presence  of 
vertical  layers  (e.g.,  ground  flora,  shrub,  high  ceuiopy)  is  recorded.  These  vertical 
layers  comprise  the  third  level  of  habitat  diversity.  A  fourth  layer  of  qualifiers 
(e.g.,  conifer,  deciduous)  then  describes  them  further. 

Different  situations  and  studies  may  have  different  habitat  classification 
schemes.  Once  the  structure  and  habitat  diversity  is  defined,  the  next  question 
is  how  to  measure  the  structure  and  habitat  diversity  Methods  for  measuring 
species  diversity  can  also  be  used  for  measuring  structm*e  and  habitat  diversity 
(Magurran  1988).  A  rather  different  approach,  differentiation  diversity,  is  re¬ 
quired  when  we  wish  to  ascertain  how  species  numbers  and  identifiers  differ  be¬ 
tween  communities  or  along  gradients  (Magurran  1988). 

The  true  value  of  studying  diversity  is  in  its  application.  It  is  believed  that  di¬ 
versity  is  a  good  indicator  of  the  well-being  of  an  ecosystem  (Magurran  1988). 
Diversity  measures  have  potential  applications  to  two  main  ideas.  First,  is  the 
idea  of  natural  resources  conservation,  which  is  underpinned  by  the  notion  that 
species-rich  commxinities  are  better  than  species-poor  ones.  Second  is  the  idea  of 
environmental  monitoring,  which  assumes  that  the  adverse  effects  of  pollution  or 
disturbance  will  be  reflected  in  a  reduction  in  diversity  (or  by  a  change  in  the 
shape  of  the  species  abundance  distribution  [Magurran  1988]).  Many  research¬ 
ers  (Bechtel  and  Copeland  1970;  Schafer  1973;  Rose  1978;  Gray  and  Mirza  1979; 
Yapp  1979;  Wu  1982;  Usher  1986;  Tomascik  and  Sander  1987)  have  shown  suc¬ 
cessful  applications  of  diversity  measures. 

It  shoidd  be  noted,  however,  that  all  of  the  previously  mentioned  diversity  meas¬ 
ures  have  some  limitations.  Species  richness  gives  only  the  nximber  of  species  or 
species  density.  Diversity  indices  based  on  the  proportional  abundance  of  species 
represent  the  number  of  species  and  their  relative  abimdance  as  a  single  figure. 
However,  they  lose  information  about  relative  species  abundance.  Species  abun¬ 
dance  models  give  the  fullest  description  of  species  distribution.  Once  the  spe¬ 
cies  distribution  of  a  commxmity  is  determined,  diversity  indices  can  be  calcu- 
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lated  from  the  distribution  of  species.  However,  species  abundance  models  also 
have  problems  in  application. 

Although  there  are  many  mathematical  distribution  models  related  to  species 
abundance,  two  models  have  foctxsed  on  the  limited  distributions  of  species:  the 
log-normal  and  the  log-series  models.  The  log-normal  model  is  said  to  fit  a  wide 
variety  of  communities  that  are  stable  and  contain  a  large  number  of  species, 
while  the  log-series  model  is  said  to  be  applicable  to  small,  stressed,  or  tmstable 
communities.  Although  the  log-normal  and  the  log-series  models  have  been 
widely  used  and  have  good  theoretical  explanations,  there  are  three  problems 
with  using  these  models.  First,  the  majority  of  real  commiinities  may  not  be 
satisfactorily  described  by  either  the  log-normal  or  the  log-series.  It  is  very  rare 
to  find  a  commiinity  with  the  exact  log-normal  distribution  or  log-series  distribu¬ 
tion.  Second,  either  the  log-normal  or  the  log-series  distribution  may  result  as 
an  artifact  of  the  sampling  procedure.  Combining  smaller,  separate  samples 
may  produce  a  log-normal  distribution,  even  if  the  individual  samples  each  show 
a  log-series  distribution.  On  the  other  hand,  a  small  sample  taken  from  a  log¬ 
normal  community  may  produce  a  distribution  resemhhng  a  log-series.  Third, 
even  if  the  log-normal  or  the  log-series  distributions  fit  the  community  very  well 
by  the  ‘goodness  of  fit’  test,  the  fitted  log-normal  or  the  log-series  model  does  not 
provide  a  good  estimation  of  the  abundance  of  the  most  prevalent  species  or  the 
least  prevalent  species.  This  is  because  the  most  abimdant  and  least  abimdant 
species  are  represented  at  the  tails  of  the  distribution.  Most  biodiversity  meas¬ 
ures  like  the  Shannon  index  and  the  Simpson  index  depend  on  a  large  niimber  of 
rare  species.  Also,  the  quality  of  a  community  is  mainly  determined  by  the  most 
abimdant  species.  Therefore,  using  a  theoretic  distribution  in  biodiversity 
analysis  causes  a  precision  problem. 

Although  the  problems  discussed  here  are  critical  in  a  diversity  study,  they  have 
not  been  solved  yet.  This  study  intends  to  explore  a  new  method  of  measuring 
diversity  and  fill  the  gap  in  the  recent  diversity  study.  A  new  stochastic  dynam¬ 
ics  model  will  be  developed  to  model  both  the  deterministic  population  growth 
and  population  fluctuations.  This  model  will  use  a  variable  number  of  species 
and  species  abundance  to  calculate  any  kind  of  diversity  indices.  Because  both 
the  dynamic  changes  and  stochastic  fluctuations  are  included  in  the  new  model, 
the  stochaustic  dynamics  model  has  the  potential  for  extensive  use  in  natural  re¬ 
sources  management  and  environmental  monitoring.  Future  applications  of 
these  stochastic  dynamics  models  include  (1)  providing  standard  diversity  meas¬ 
ures,  (2)  monitoring  the  development  of  plant  communities  in  terms  of  species 
diversity  and  structure  diversity,  (3)  testing  the  significance  of  the  influence  of 
human  activities  on  plant  communities,  and  (4)  estimating  rehabilitation  time 
for  a  disturbed  plant  community. 
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3  Stochastic  Dynamic  Models 

Although  the  real  mechanisms  for  the  development  of  a  plant  community  are 
very  complicated,  mathematical  models  can  replace  the  complex  biological  reality 
with  some  idealized  hypothetical  systems.  In  fact,  many  simphfied  mathemati¬ 
cal  models  can  interpret  complex  consequences  and  predict  the  behavior  of  elabo¬ 
rate  natxiral  systems.  In  this  report,  we  will  discuss  several  population  dynam¬ 
ics  models  with  environmental  stochasticity  and  demographic  stochasticity. 

A  stochastic  dynamics  model  can  be  derived  by  manipulating  a  deterministic 
model  to  include  demographic  stochasticity  and  environmental  stochasticity.  A 
deterministic  model  describes  the  mechanisms  controUing  the  population  growth 
and  decay.  The  values  of  the  parameters  in  the  deterministic  model  are  assumed 
to  be  known.  Whenever  the  initial  conditions  are  given,  the  deterministic  models 
provide  exact  predictions  of  future  populations.  In  the  real  world,  however,  few 
communities  have  the  population  dynamics  described  by  deterministic  models. 
Instead,  populations  of  most  natmal  communities  have  fluctuating  growth  be¬ 
cause  of  demographic  and  environmental  stochasticity  (or  noise). 

Turelli  (1986)  provides  definitions  of  environmental  stochasticity  and  demo¬ 
graphic  stochasticity.  Demographic  stochasticity  (or  within-individual  variabil¬ 
ity)  is  the  variation  of  individuals  who  appear  to  be  identical  but  have  different 
life  lengths  and  produce  different  numbers  of  offspring.  Integer-valued  stochas¬ 
tic  models  are  t3rpically  used  to  investigate  the  consequences  of  the  demographic 
stochasticity.  Environmental  factors  vary  tmpredictably  through  time  in  ways 
that  affect  all  individuals.  This  variation  is  called  environmental  stochasticity. 
Most  analyses  of  the  consequences  of  environmental  stochasticity  begin  by  add¬ 
ing  a  noise  term  in  the  deterministic  model  (May  1973;  CapocelH  and  Ricciardi 
1974;  TuckweU  1974;  Goel  and  Richter-Dyn  1974;  Turelli  1977).  This  produces 
stochastic  difference  and  differential  equations  with  continuous  ranges. 

In  this  paper,  we  first  introduce  a  d5Tiamics  model  of  population  growth.  This 
model  derives  birth  and  death  rates  as  they  relate  to  population  growth  from  the 
relationships  among  plants.  It  also  derives  these  rates  from  the  relationship 
between  plants  and  the  environment.  The  model  will  serve  as  the  deterministic 
part  of  the  stochastic  model.  Next,  we  present  some  well-known  stochastic 
models  such  as  birth  and  death  processes  and  a  diffusion  process.  In  this  part, 
we  demonstrate  different  methods  to  solve  a  variety  of  stochastic  differential 
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equations.  Next,  we  will  show  how  to  incorporate  both  demographic  noise  and 
environmental  noise  into  a  single  model  that  considers  the  joint  elfects  of 
demographic  and  environmental  stochasticity.  Finally,  we  will  apply  the  model 
to  predict  the  dynamics  of  the  understory  in  a  bottomland  forest. 


A  Simulation  Model  of  Population  Dynamics 

The  dynamics  of  the  population  growth  of  a  community  can  be  mathematically 
described  if  the  functional  behavior  of  the  rate  of  growth  is  known.  The  litera¬ 
ture  contains  many  different  representations  of  the  growth  rate  (Nisbet  and 
Gurney  1982;  Streifer  1974;  Hallam  1986).  Most  of  these  representations  have 
two  basic  factors  in  common:  survivorship  and  recruitment.  Survivorship  can  be 
expressed  by  the  percentage  of  individuals  that  survive  from  one  time  period  to 
the  next.  Recruitment  is  the  addition  of  new  individuals  to  the  population  by 
immigration.^  and  births.  For  plant  communities,  the  emigration  can  be  ne¬ 
glected  because  plants  usually  do  not  emigrate  once  they  become  estabhshed  in 
the  community. 

A  Discrete  Model  of  Plant  Population  Dynamics 

Consider  a  1-hectare  plant  community  consisting  of  10,000  space  units  (i.e.,  1 
space  unit  =  1  m^).  Each  individual  of  i-th  species  occupies  A‘  space  unit.  Thus, 
the  space  occupied  by  all  plants  at  time  t  is  =  'Lxl*A!,  where  is  the  niimber 

of  individuals  of  the  i-th  species  at  time  t. 

The  change  in  the  number  of  individuals  of  the  i-th  species  from  time  t  to  time 
t-i-1  depends  on  the  number  present  at  time  t,  the  nximber  recruited  (immigration 
and  birth),  and  the  nmnber  of  individuals  that  survive.  This  relationship  is  ex¬ 
pressed  as: 

(D*  =  s‘  (xi  +  F) 

where  x^^^  £md  x^  are  the  number  of  individuals  of  the  i-th  species  at  time  t-nl 
and  t,  respectively,  s'  is  the  net  survivorship,  and  r'  is  the  net  recruitment. 


Equations  in  each  chapter  begin  with  (1). 
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Competition,  insects  and  disease,  senescence,  and  occasional  catastrophe  (natu¬ 
ral  catastrophe  and  human  disturbance)  are  the  most  important  sources  of  mor¬ 
tality  for  terrestrial  plants.  Mortality  (or  survivorship)  is  incorporated  into  the 
model  in  three  ways.  First,  the  inter-specific  siirvive  rate  (S)  is  included  in  the 
model  as  a  given  number  ranging  from  0  to  1.  The  inter-specific  survival  rate  of 
the  i-th  species  (S')  represents  the  intrinsic  survivability  of  the  i-th  species.  Sec¬ 
ond,  to  account  for  the  mortality  due  to  a  catastrophe,  the  inter-specific  survival 
rate  is  multiplied  by  the  catastrophe  index  C  (0<C<1).  We  assume  that  all  spe¬ 
cies  suffer  from  the  catastrophe  but  their  tolerances  are  different.  Third,  the 
risk  of  mortality  from  competitive  interference  is  likely  to  increase  with  an  in¬ 
crease  in  the  community  density.  In  this  model  we  use  a  discrete-time  analog  of 
the  Lotka-Volterra  competition  equation  which  is  expressed  as: 

(2)  Qi  =  gi(5-^) 


where  K  is  the  maximvun  nvunber  of  individuals  of  the  i-th  species  emd  ajj  are 

the  competition  coefficients.  We  assume  that  for  each  i,  g(x)>0,  and  g(x)'<0  for 
x>0  (i.e.,  g(x)  is  a  positive  decreasing  function).  Thus,  the  survivorship  compo¬ 
nent  of  equation  (1)  is  now  written  as: 


(3) 


S' *C 
l  +  Q‘ 


The  number  of  individuals  recruited  depends  on  the  inter-specific  recruitment 
potential  (R'),  the  availability  of  space  (A/A)  at  time  t  and  the  immigration  (I). 
The  inter-specific  recruitment  rate  represents  the  rate  of  addition  of  population 
by  birth.  As  the  community  develops  and  as  space  becomes  occupied,  the  poten¬ 
tial  for  further  recruitment  is  proportionally  reduced.  Therefore,  the  inter¬ 
specific  recruitment  potential  (R')  should  be  multipHed  by  a  space  constraint  (V) 
which  is  a  function  of  the  ratio  of  space  occupied  and  total  space,  i.e.,  V  =  f(A/A). 
The  combined  recruitment  rate  is: 

(4)  r‘  =  (R‘:c;  +  r)V 

Putting  s'  and  r'  in  equation  (1),  we  obtain  the  final  version  of  the  dynamics 
model. 

<5)  ^Li  =  +  r)V] 

1+  y 
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This  model  produces  a  “community”  with  a  variable  number  of  species  and  vari¬ 
able  population,  which  can  then  be  used  to  calculate  diversity  indices. 


The  Continuous  Model  of  Plant  Population  Dynamics 


When  the  population  size  of  a  commxmity  is  large,  a  continuous  model  can  ap¬ 
proximate  population  growth.  Let  and  x‘  respectively  be  the  numbers  of  i- 

th  species  at  times  t-i-dt  and  t.  The  change  in  the  number  of  individuals  of  the  i- 
th  species  from  time  t  to  time  t+dt  is: 


qi  qi  qi  *  c 

(6)  4,,  -  x;  =  {xiiB!  V  -  (1-?-^)]  +  r  v}*dt 


1+  Q* 


1+Q‘  1+Q' 


=  [x‘(fr-d‘)-Hm‘]*dt 


Dividing  dt  on  both  sides  of  (6)  and  letting  dt  go  to  0,  we  have  the  following  dif¬ 
ferential  equation. 


(7) 


—(xi)  =  xi(h'-d)+  m‘ 


=  jc'*h‘  -1-  m’ 

*  c 

where  b‘  =(R' - r  V)  is  the  birth  rate,  d‘  =(1-- — )  is  the  death  rate,  m‘ 

1  +  Q'  1+Q‘ 

is  the  net  immigration  rate  and  h'  is  the  net  growth  rate.  The  solution  of  equa¬ 
tion  (7)  is  an  exponential  growth  function  with  immigration. 

(8)  x‘  =  (Xq  -t-  — )*exp(h’t)  -  — 

r  r 

When  h>0,  immigration  causes  sm  exponential  increase  in  popiilation  growth. 
When  h<0,  immigration  causes  an  exponential  decrease  in  poptJation  growth, 
and  there  is  high  competition  in  the  community.  The  new  immigrants  increase 
the  competition  and  speed  up  the  decrease  of  popiilation  growth.  As  in  equation 
(5),  equation  (8)  also  produces  a  “community”  with  a  variable  number  of  species 
and  variable  population. 
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Demographic  Stochasticity 

As  discussed  in  the  dynamics  models,  the  birth  rate,  death  rate,  and  immigration 
rate  are  assumed  to  be  constant.  In  the  natvu-al  world,  however,  this  assumption 
does  not  usually  hold  true.  Individuals  of  identical  plant  t3T)es  may  have  differ¬ 
ent  life  lengths  and  produce  different  numbers  of  offspring.  This  variation 
among  the  individuals  is  called  demographic  stochasticity.  There  is  vast  litera¬ 
ture  on  the  modeling  of  demographic  noise.  TureUi  (1986)  gives  a  very  good 
overview  of  demographic  stochastic  models. 

In  this  paper  we  will  use  the  birth  and  death  processes  to  describe  the  demo¬ 
graphic  stochasticity.  The  birth-death  processes  describe  population  dynamics 
with  biologically  accurate,  interpretable  birth  and  death  rates  and  are  applicable 
to  individual  numbers  of  every  size.  First,  we  start  with  an  ‘external  birth’ 
model,  the  Poisson  process  of  immigration. 

The  Poisson  Model  of  Immigration 

Suppose  the  chance  of  an  event  (immigration)  in  the  small  time  interval  (t,  t+dt) 
is  >,dt+o(dt),  where  the  last  term  is  a  remainder  term  which  becomes  negligible 
compared  with  dt  as  dt  gets  smaller,  and  may  be  consequently  neglected  in  com¬ 
parison  with  the  first  term.  This  chance  Xdt  is  assumed  to  be  independent  of  the 
number  of  previous  happenings,  and,  moreover,  each  event  is  assumed  to  be  in¬ 
dependent.  Therefore  the  chance  of  two  events  (or  more)  occurring  in  the  time 
interval  (t,  t-i-dt)  is  o(dt),  and  is  also  negligible.  Let  the  probability  of  r  events 
(and  no  more  nor  less)  in  the  time  interval  (0,t)  be  p^.  Then  we  can  mathemati¬ 
cally  represent  the  whole  set  or  distribution  of  probabilities  (r=0,  1,  2,  ...)  by 
the  generating  function  in  0, 

(9)  n(0)  =  i  =  0  Pr0i 

As  the  Pr  are  dependent  on  t,  we  note  that  11(0)  is  also  a  function  of  t,  and  shall 
write  it  more  fiilly  as  n/0).  Then  by  the  ndes  of  probability,  as  the  increase  in 
the  total  number  of  events,  during  the  further  interval  dt,  is  independent  of  the 
previous  total  number  occurring  in  the  interval  (0,  t),  it  follows  that 

=  n.(0){l-Xdt-f)i0dt}  or 

(10)  J-log(/^(9)=M0-l) 


whence,  as  IlofO)  =  1, 
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(11)  n/G)  =  exp[Xt(0-l)],  and 

(12)  Pr  =  (^t)r*exp(-X.t)/r! 

Equations  (11)  and  (12)  specify  the  Poisson  distribution  with  mean  m=A,t.  When 
X  is  time-varying,  the  mean  is  obtained  by  the  following  equation: 

t 

(13)  m  =  J  X(u)du 

0 

The  Simple  Birth  and  Death  Process 

We  assume  that  the  probability  for  any  given  individual  to  give  birth  in  the  time 
interval  (t,  t-f-dt)  is  Idt  while  the  probability  of  dying  in  that  time  interval  is  mdt. 
Equivalently,  denoting  by  dX(t,  t-^dt)  the  increment  X(t-t-dt)-X(t)  of  the  population 
size  in  (t,  t-ndt)  we  can  make  the  following  assmnptions: 

(14a)  P{dX(t,  t-i-dt)  =  1/X(t)  =  n}  =  Indt  +  o(dt) 

(14b)  P{dX(t,  t-i-dt)  =-l/X(t)  =  n}  =  mndt  +  o(dt) 

(14c)  P{  I  dX(t,  t-i-dt)  I  >  1  /  X(t)  =  n}  =  o(dt) 

From  (14)  it  follows 

(15)  p{dX(t,  t+dt)  =  0/X(t)  =  n}  =  l-(l-fm)ndt  +  o(dt) 

From  (14)  and  (15)  we  easily  obtain: 

p^(t+dt)  =  p^(t)  [l-(l-i-m)ndt]  -i-  p^.,(t)(n-l)ldt  -i-  p__^j(t)m(n-i-l)dt  +  o(dt) 
or,  the  hmit  as  dt  approaches  0: 

(16)  -  -n(l-fm)p„(t)  -1-  l(n-l)p„.i  (t)  +  m(n+l)p^^j(t) 


with  the  initial  condition: 


P„(0)= 


hn=j 

0,n^j 


(17) 
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Before  showing  how  equation  (16)  can  be  solved  to  determine  the  functions  p„(t), 
let  us  calculate  the  mean  population  size  E[X(t)]  and  its  variance  Var  {X(t)}  by  a 
straightforward  procedure.  Use  the  definition: 

(18)  E[X(t)]  =  n  =  0  np„(t) 

and  differentiate  both  sides  with  respect  to  t: 

(19)  =  n  =  0  n4-p„(t) 

dt  dt 

=-n  =  0  n2(A,-hp)p„(t)-i-Xn  =  0  n(n-l)p„.(t)H-pn  =  0  n(n+l)p^^j(t) 

where  the  last  equality  follows  from  (16).  Setting  n-l=n  and  n4-l=n  in  the  second 
and  third  sums  on  the  right-hand  side  (r.h.s.)  of  (19),  respectively,  after  some 
straightforward  algebra  we  are  led  to  the  following  differential  equation: 

(20)  =  (^+li)E[X(t)] 

On  the  other  hand  from  (19)  and  (20)  we  obtain  the  initial  condition 

(21)  E[X(0)]=j 

From  (20)  and  (21)  we  immediately  get 

(22)  E[X(t)  /  X(0)  =j]  =j*exp[(X-p)t] 

By  a  similar  procedure  an  equation  for 

(23)  Var  {X(t)}  =  E[X2(t)]  -  {E[X(t)]}2 

can  be  derived  and  solved  with  the  initial  condition 


Var  {X(0)}  =  0 


We  have  the  following  result: 


(24)  Var{X(t)/X(0)=j}  =  ^ 


X=  n 


j{X+  fj) 


[  X- ^ 
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Note  that  in  the  case  this  variance  not  only  depends  on  the  intrinsic  “growth 
rate”  X-p  but  also  on  the  sum  X+p.  The  sum  X+p=R  defines  the  noise  strength, 
called  demographic  noise  (Nisbet  and  Gurney  1982). 

Now  we  turn  to  the  solution  of  equation  (16).  We  first  introduce  the  monument 
generating  function,  M(0,  t),  defined  as  follows  in  terms  of  the  p„s: 

(25)  M(0,  t)  =  n  =  0  e®“p„(t) 


where  p  is  a  dummy  variable.  Multiplying  both  sides  of  equation  (16)  by  e®“  and 
summing  over  n  fi'om  0  to  oo,  one  can  easily  obtain: 


<26) 

dt  oO 


with  the  initial  condition 


(27)  M(0,  0)  =  e®' 


The  general  solution  of  equation  (26)  can  be  found  by  the  method  of  characteris¬ 
tics. 


(28) 


Xv{e,t)-\_ 

l-(Ar-l){e^-l)T 


/d 

,X  =  n 


where 

Now  we  define  function  F(s,  t)  as  the  probability  generating  function: 


(29)  F(s,  t)  =  n  =  0  p„(t)s“  =  MGog(s),  t) 


(30) 


F{s,t) 


//(l-Q:)-(A-//g)iy  ’ 
H-Xa-X{\-a)s 
'i-(;h-iX5-i)T 

i-;u(5-i) 


,x=n 


where  a=e“’'^’‘ 
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The  meaning  of  ^-|a  is  intuitive.  It  is  the  net  rate  of  population  change.  To  un¬ 
derstand  the  meaning  of  A.+p  we  introduce  the  diffusion  approximation.  We  re¬ 
write  equation  (16)  as 

(31)  =  -[X(n)+p(n)]p_^(t)-f)i(n-l)p  (t)-i-p(n-i-l)p„^j(t) 

dt  "■i 

We  use  Taylor  expansions  to  approximate  the  various  functions  of  n±l  in  (31). 
This  involves  treating  n  as  a  continuous  variable  and  reinterpreting  p^(t)  as  a 
probability  density.  Discarding  all  terms  in  the  Taylor  expansions  that  are  of 
third  or  higher  order,  we  have  the  approximate  equation  of  (31), 

(32)  =  -ir[(A(n)  (<)] +m)P.  (<)] 

With  (32)  it  is  easy  to  show  that 

(33)  E(dn)  =  [?i(n)  -  p(n)]dt 

(34)  E[(dn)^]  =  [Mn)-Kp(n)]dt 

Equation  (34)  explains  why  we  call  X-t-p  the  strength  of  demographic  noise. 

Now  we  have  established  simple  closed  form  expressions  for  the  probability  p„(t) 
and  calculated  explicitly  mean  and  variance  of  the  population’s  size.  All  this  al¬ 
lows  us  to  draw  an  accurate  picture  of  the  population’s  time  evolution.  In  par¬ 
ticular,  we  can  easily  get  information  about  the  extinction  probability  of  the 
population  by  taking  the  limit  of  (30)  as  s  goes  to  0. 

(35)  =  f  , 

Xi 

_1  -t-  Xt 

Environmental  Stochasticity 

Environmental  fluctuations  inevitably  produce  fluctuations  in  popiilation  levels. 
A  general  question  is  how  species  dynamics  and  interactions  translate 
environmental  fluctuations  into  temporal  and  spatial  patterns  of  popiilation 
abundance.  Because  there  is  a  lack  of  mathematical  machinery  available  to 


,A  =  // 


CERL  TR  00-5 


26 


analyze  nonlinear  multidimensional  stochastic  processes,  environmental 
fluctuations  are  incorporated  into  a  deterministic  model  in  a  way  that  is 
biologically  meaningful  yet  mathematically  tractable.  The  most  common  way 
used  in  the  literature  is  to  add  a  noise  term  in  a  differential  or  difference 
equation  (May  1973;  Capocelli  and  Ricciardi  1974;  Tuckwell  1974;  Goel  and 

Richter-Dyn  1974;  Turelh  1977).  Here  we  discuss  the  two  most  common  types  of 
stochastic  differential  equations:  =  h(x)  +  k(x)A(t)  and  dx  =  h(x)dt  +  k(x)dW. 

Stochastic  Differential  Equations  (1) 

f  =  h(x)-fk(x)A(t) 


We  first  deduce  the  so-called  kinetic  equation  from  Markovian  property.  The  ki¬ 
netic  equation  is  the  general  form  of  many  stochastic  differential  equations. 

Consider  a  Markov  process  with  a  continuum  of  state  values  in  continuous  time. 
Its  transition  probability  density  function  (p.d.f.),  I  (x,  t  I  xq,  tq),  satisfies  the 

Chapman-Kohnogrov  equation  (Bartlett  1966;  Ross  1983;  Bharucha-Reid  1960): 

(36)  fix,  1 1  Xo,  t„)  =  Jdy/(x,  1 1  y,  T)/(y,  x  I  x,,,  t^) 

with  t>z>tg  arbitrary  instants  and  X(t)=x,  X(T)=y,  X(t5)=x„.  Equation  (36)  is  to  be 
looked  at  as  a  compatibility  relation  holding  for  any  Markov  process,  but  it  is  not 
sufficient  to  determine  the  process’  transition  p.d.f.  To  accomplish  this  task,  fur¬ 
ther  assumptions  besides  the  Markov  assiunption  are  necessary.  First  let  us  re¬ 
write  equation  (36)  in  a  differential  form, 

(37)  fix,  t+At  I  x^,  t,)-/(x,  1 1  x„,  to)  =  Jdy/(x,  t-^At  I  y,  t)/(y,  1 1  x^,  to)-/(x,  1 1  x^,  t^,) 


Let  us  now  consider  an  arbitrary  function  R(x)  vanishing  at  the  end  points  of  the 
state  space  sufficiently  rapidly,  together  with  its  derivatives  of  all  orders.  Multi¬ 
plying  both  sides  of  (37)  by  R(x)/At  and  integrating  over  the  state  space,  we  ob¬ 
tain: 


Jcic/?(x) 

(38)  =ij’ cb:R(x)Ji/y/’(x,/  +  Atly,t)f(y,t|xo,fo) 
— chR(x) /’(x,dxo,/o) 


Substituting  the  Taylor  expansion  about  the  point  y  for  R(x)  in  the  first  integral 
on  the  right-hand  side  of  equation  (38): 
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(39) 


R(x)  =  R(y)+'E^ 


Rjy)  jx-y)" 


and  taking  the  limit  as  At  goes  to  0,  we  obtain: 

|jxR(x)|-=  lim-^  JdyR(y)f(y,t|Xo,to)  ^xf(x,t  + At|y,t) 

At->0 

+  iy{^f(y.t|xo,to) 


(40) 


n=i 


X  lim  fclx(x  -  y)" f(x,t  +  At|y,t)} 

At->-0 

-lim-sf  JdxR(x)f(x.t|Xo.to) 


or: 


(41) 


]dxR{x)-^ 


^  fix,t\xo,to)A„ix,t) 


n-] 


having  used  the  normalization  condition: 

(42)  \dx:f  {x,t  + ixt\y,t)  =  \ 
and  having  set: 

(43)  A„(x,t)  =  \lV[i^\^yiy  ~  x)"  f{y,t  +  A/|x,/)>  (n=l,  2, ...) 

The  integration  by  parts  of  the  r.h.s.  of  equation  (41),  in  which  the  vanishing  of 
R(x)  and  its  derivatives  at  the  ends  of  integration  interval  is  used,  shows  that: 


(44) 


J  f  {y,t\xQ,h)A„{y,t) 

=  (-!)"  ]^xR{x)i,[Mx,t)nxAxoM 


Equation  (41)  thus  yields: 


(45)  Jcici?(x){-^-|;^^[^(x,/)A^,^l^o,fo)]}  =  0 
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Now  due  to  the  arbitrariness  of  the  function  R(x),  the  bracketed  terms  must  be 
identically  zero,  and  we  have  our  desired  result: 
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(46) 

n=\ 

This  is  the  kinetic  equation  holding  under  the  sole  assumption  that  the  process 
under  consideration  is  Markovian.  The  functions  A_j(x,  t)  defined  by  (46)  are 
csdled  infinitesimal  moments  of  the  process. 

Consider  a  poptdation  growth  model  with  a  stochastic  growth  rate  due  to  the  en¬ 
vironmental  fluctuations: 

(47)  r  =  “r  +  A(t) 


where  A(t)  is  a  noise  term  due  to  the  environmental  stochasticity  and  ”r  is  a  de¬ 
terministic  net  growth  rate.  Then  we  have  a  simple  stochastic  differential  equa¬ 
tion: 

(48)  —  =  x*”r  +x  *A(t) 

dt 


Extending  the  situation  above,  we  consider  a  general  linear  equation  of  the  type 

(49)  —  =  h(x)  +  k(x)A(t) 

dt 


where  h  and  k  are  assigned  functions  and  A(t)  is  a  stochastic  process.  Clearly, 
the  solution  of  (49),  x(t),  is  a  random  function.  Its  determination  cannot  be  ac¬ 
complished  unless  A(t)  is  specified. 

Let  us  assume  that  A(t)  in  (49)  is  a  stationary  process  with  a  0  mean  and  with  a 
rather  narrow  and  peaked  correlation  function: 

(50a)  E[A(t)]  =  gl  =  0 

(50b)  E[A(ti)A(t2)]  =  g2(tj,  t^)  =  g2(t2-ti) 

where  g^ft)  is  appreciably  non-zero  only  in  the  neighborhood  of  t=0  with  a  very 
sharp  maximum  at  t=0.  More  generally,  for  any  group  of  instants  t^,  t^,  ...,  t^  all 
lying  close  to  each  other  we  set: 


(51) 


E[A(t,)A(t2) ...  A(t„)]  =  g„(tj,  t„) 
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and,  again,  assume  that  the  n-th  order  correlation  function  g„  has  a  sharp  maxi¬ 
mum  at  tj=t2=...=t„,  being  otherwise  effectively  0.  Finedly,  we  assume  that  when 
tj,  tj,  ...,  t^  are  proximal  to  each  other,  and  also  when  t^+1,  t+2, ...,  t^  are  proximal 
but  far  from  the  group  tj,  t^  and  so  on,  then: 

(52)  E[A(t,)  ...  A(t>(t +1) ...  A(t>(t +1) ...  A(t^) ...  ] 

=  E[A(t,) ...  A(t,)]E[A(t,-t-l) ...  A(t,)]E[A(t,-i-l) ...  A(g] ... 

=  -  t,)gA+l  t,)gp(t,-Hl ...  V  ... 

where  the  functions  g^  have  already  been  qualitatively  specified. 


All  these  assumptions  about  the  stochastic  process  Aft)  appearing  in  equation 
(49)  may  look  rather  artificial  at  this  stage,  but  the  motivation  for  them  will  soon 
be  apparent.  With  this  in  mind,  let  us  perform  a  change  of  variable  in  (49)  by 
setting: 

(53)  y  =  0(x)  x=a)-l(Y) 
with 

(54)  <I>(x)=}^ 


Then,  equation  (49)  is  changed  into: 


(55) 


dy 

dt 


H{y) 

Kiy) 


+  A(t) 


upon  setting: 

(56)  H(y)  =  h[a)-1(Y)] 

K(y)  =  k[cb-l(Y)] 

The  advantage  of  this  procedxme  is  that  we  have  constructed  a  stochastic  process 
y(t)  defined  by  the  simpler  equation  (55)  in  which  Aft)  appears  in  a  pirrely  addi¬ 
tive  way.  Due  to  the  above  assumptions  on  Aft),  we  now  expect  y(t),  and  hence 
x(t),  to  be  Markovian.  Its  transition  p.d.f.  f^fy,  t/y^)  thus  satisfies  the  kinetic 
equation  which  we  have  derived  in  (46). 
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(57) 


dfy 

dt  ^  n\  dy" 


To  evaluate  the  infinitesimal  moments  we  first  express  the  increment  of  y 
over  a  small  time  interval  dt  by  means  of  the  approximation: 

H(v)  '7' 

(58)  dy  =  y(t-Kdt)  -  y(t)  «  _;  -:  dt  +  [  A(T)dT 

K(y) 


where,  here  and  in  the  following,  the  value  of  the  process  at  time  t  is  considered 
as  fixed.  Note  that  equation  (58)  requires  that  H(y)  and  K(y)  be  smooth;  how¬ 
ever,  the  smoothness  of  the  sample  paths  of  A(t)  is  not  implied.  Taking  the  ex¬ 
pectation  of  both  sides  of  (58),  due  to  (50),  in  the  limit  as  dt  0  we  obtain: 


(59) 


df->0  01 


Hiy) 

Kiy) 


To  calculate  we  now  square  both  sides  of  (58)  and  obtain: 


(60) 


jjf  \  t-^dt  t+dt  t+d: 

(dy)2 « 0[(dt)2]  +  2dt^^^  f  A(T)dT-i-  [  A(T)dT  f  A(0)de 

K(y)  •'  •'  •' 


Upon  taking  the  expectations  and  after  dividing  by  dt,  for  small  dt  we  are  left 
with: 

t+dl  t+dt 

(61)  B,dt«E[(dy)1=  j  Adt  J  g,(T-0)d0 

t  t 

after  making  use  of  equation  (50).  Using  the  earlier  specified  qualitative  behav¬ 
ior  of  gj,  it  then  follows: 

(62)  Bj  « <3  with  =  J” 

and  with  the  result  becoming  exact  in  the  limit  as  dt  — >  0. 

Proceeding  along  similar  fines,  it  is  not  difficult  to  become  convinced  that  due  to 
the  assumed  properties  (50)  and  (52)  of  A(t)  the  following  relationship  holds: 


(63) 


dt*B„(y) «  E[(dy)”]  =  o(dt)  (n  =  3,  4, ...) 
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Equation  (57)  thus  becomes  a  forward  diffusion  equation: 


(64) 


Idy'^ 


with  Bj(y)  and  a*  given  by  (59)  and  (62),  respectively.  The  conclusion  is  that 
equation  (55)  can  be  thought  of  as  defining  a  diffusion  process  y(t)  whose  drift 
equals  the  deterministic  part  of  the  r.h.s.  of  the  equation  while  its  infinitesimal 
variance  depends  exclusively  on  the  characteristics  (gjfi:))  of  the  random  part  of 
the  equation.  Fiirthermore,  if  we  impose  that  p{y(0)=yj=l,  with  y^,  being  non- 
random,  then  the  specification  of  such  a  diffusion  process  is  rmique. 


Let  us  now  examine  the  infinitesimal  moments  A^^fx)  of  the  Markov  process  x(t) 
defined  by  equation  (49).  Denoting  its  transition  p.d.f  by  f^(x,  t/x„),  we  know  that 
for  small  dt  we  have: 


(65)  A^(x)dt « J(x'-x)'’f^(x',  dt/x)dx’ 

=  [a)'(y’),  dt/x]d[o'(y)] 

X 

=  l[d>  ‘(y')-x]''f^[y’,  dt/0(x)]dy'  (n  =  I,  2, ...) 


having  made  use  of  the  one-to-one  transformation  (53)  between  the  transition 
p.d.f.s  of  the  processes  x(t)  and  y(t).  Let  us  now  expand  0  *(y')  as  a  Taylor  series 
about  the  point  y=d)(x): 


(66) 


d>  '(y')  =  C>'^[(I)(x)]  -I-  aj(x)[y'-<l)(x)]  -i-  2  a2(x)[y'-(I>(x)]^ 
o  I 

n  =  3  — — ,n! 


with  a„(x)  = 


d’'F-\y' 


idy  ly=<j)(x) 


It  is  easy  to  see: 

(67)  ax=k(x) 


a2=k'  (x)k(x) 
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where  k'  (x)=dk(x)/dx 

Using  (66)  and  (67)  and  taking  the  hmit  dt  0,  we  obtain: 

1  a  (x^ 

(68)  A,(x)  =  k(x)B,[0(x)]  -h  -k’(x)k(x)B,[<D(x)]  -h  Z  B  J(D(x)] 

2  «=3  n\ 

Making  use  of  equations  (59)  and  (61)  we  thus  find: 

<>2  ^fc2/ 

(69)  Ai(x)  =  h(x)  +  — - ,dx 

Using  the  same  procedure,  we  can  have: 

(70)  A  (x)  =  a^k^(x) 

2 

A„(x)  =  0  (n=3,  4, ...) 

Thus,  we  conclude  that  the  Markov  process  x(t)  defined  by  equation  (49)  is  a  dif¬ 
fusion  process  with  drift  and  infinitesimal  variance  given  by: 

dk^ix 

(71)  ApO  =  h(x)  +  — - ,dx 


Ap)  =  a^^(x) 

Stochastic  Differential  Equations  (2) 
dx  =  h(x)dt  +  k(x)dW 

In  plant  communities  the  effect  of  environmental  fluctuations  is  cumulative  with 
plants  through  time.  Although  the  environmental  stochastic  influence  at  each 
time  point  can  be  considered  as  some  noise,  the  cumtdated  realistic  noise  could 
be  some  other  process.  This  case  can  be  described  by  the  following  equation: 

(72)  dx  =  h(x)dt  +  k(x)dW 

where  W  is  some  stochastic  process  such  as  Brownian  motion  or  Wiener  process. 
Equation  (72)  is  called  an  Ito  equation.  This  equation  can  be  handled  by  means 
of  Ito’s  calciilus  which  differs  in  several  fundamental  ways  fi'om  classic  calculus 
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(Rozovskii  1990).  Here  we  give  one  of  the  major  results  of  Ito’s  calculus:  Ito’s 
formula,  which  can  be  used  to  find  the  solution  of  the  stochastic  model  (72). 

Ito’s  Formula:  Let  U:  [0,  oo)'R'^  be  a  continuous  function  with  continuous  par¬ 
tial  derivatives  on  (0,  oo)'R'^: 


(73) 


uo  = 


d  u 

n 


du 

dx, 


m  = 


i=l,  2, ...,  N 


d^U 


uji  =  - -  i=l,  2, ...,  N 

dX,BX, 


^  J 


and  let  Xj,  X,,  ...,  X„  be  stochastic  integrals  defined  by: 


(  I 

(74)  X,(t)  =  X(0)+J  9idW.+  J  9,d. 


then  Y  defined  by: 


(75)  Y(t)  =  u(t,  X,(t),  X,(t), ...,  X„(t)) 


is  also  a  stochastic  integral  and  its  stochastic  differential  is: 


(76)  dY=udt+y  udX;+-y  y  udxdx 

i=l  ^  i=l  i=l 


and  the  stochastic  integral  is: 


N  t 


(77)  Y(t)  =  Y(0)+  u,(s,X.,X,,...,X^)ds+X  u,(s,  X^,  X,, .. .,  X^)dXi  + 

0  i=l  0 


E  u,fe,X.,X„...,X„)cK,clX, 

^  i-1  j=i  0 


N  t 


N  t 


Y(t)  =  Y(0)+  u,d,+  2^  u^(p.dW.+  £  + 

i=l  0  i=l  0 


0 
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1  N  N  ^ 

^  i=l  j=l  0 

Three  simple  cases  of  equation  (75)  are  (1)  Y(t)=u(Wt),  (2)  Y(t)=u(t,  W^),  and  (3) 
Y(t)=u(t,  X(t)). 

Using  Ito’s  formula,  we  can  easily  find  the  solution  of  (72).  From  (76)  we  can  see 
that  the  solution  x(t)  is  also  a  diffusion  processes.  Its  drift  and  infinitesimal 
variance  are  given  by: 

(78)  Bj(x)  =  h(x) 

(79)  B,(x)  =  a^k^(x) 

To  make  this  discussion  more  concrete,  we  consider  the  following  homogeneous 
unrestricted  population  growth  model  with  environmental  stochasticity.  Other 
situations  such  as  a  density  dependence  model  can  be  solved  in  a  similar  way. 
Taking  the  deterministic  functions  h(x)  and  k(x)  in  (58)  to  be  rx(t)  and  ax(t)  re¬ 
spectively,  we  have  the  simple  model: 

(80)  dx(t)  =  rx(t)  +  ax(t)dBj 

where  r  is  the  intrinsic  growth  rate  and  B^  is  a  Brownian  motion.  To  solve  (66), 
let  Y(t)=u(x(t)=ln(x(t)),  where  x(t)>0.  Using  (62),  we  obtain  the  stochastic  differ¬ 
ential  of  u(x). 

2 

(81)  du(x)  =  (r - )dt  +  adB^ 

2 

Thus,  the  solution  for  the  simple  population  growth  model  (80)  is 

2 

(82)  x(t)  =  x(0)exp[(r - )dt  +  adB^ 

2  ^ 

The  expectation  of  x(t)  is  exp(rt)*E[x(0)].  If  x(0)  is  fixed,  the  expectation  is  the 
same  as  the  deterministic  model.  In  the  deterministic  situation  (i.e.,  x(t)  is  non- 
random),  x(t)-^  00  as  t->  00  whenever  r  is  positive.  However,  in  equation  (82)  x(t) 
may  ^  oo  or  0  even  if  r  is  positive.  If  r  is  greater  than  12,  x(t)->  oo  almost 
surely  as  t— >  co.  When  r  is  less  than  ax(t),  x(t)  goes  to  zero  almost  surely  as 
t— >  CO.  When  r  is  equal  to  ax(t),  the  poptdation  growth  is  then  completely  con¬ 
trolled  by  the  environmental  fluctuations. 
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Joint  Effect  of  Demographic  and  Environmental  Stochasticity 

In  the  real  world,  both  demographic  and  environmental  stochasticity  are  un¬ 
avoidable.  Given  the  difficulties  of  analyzing  each  separately,  it  should  come  as 
no  surprise  that  few  multispecies  analyses  incorporate  both.  To  show  how  to 
model  the  combined  effect  of  environmental  and  demographic  noise,  let  us  start 
with  the  birth  and  death  process. 

We  rewrite  the  birth  and  death  process  as: 

(83)  -  -n(A.+p)p__(t)  -I-  A.(n-l)p„,i  (t)  -f  p(n-i-l)p„-Hl(t) 

dt 


or 

,  =  -(>.„-i-|i„)p,(t)  (t)  -f-  p„,iP„.(t) 

at 

where  \  =nA.  and  p,,=np.  X,-t-p=R  is  the  strength  of  demographic  noise. 

We  suppose  that  the  deterministic  dynamics  are  density  dependent,  i.e.,  the  lo¬ 
gistic  growth: 

(84)  dX/dt  =  Xr(l-X/K) 

Environmental  fluctuations  are  included  by  adding  a  noise  term  g(X)dW  in  (84). 
The  corresponding  stochastic  equation  is  given  as: 

(85)  dX  =  Xr(X)dt  -h  g(X)dW 

where  W  is  a  stamdard  Wiener  process.  Different  choices  for  the  function  g(X) 
can  be  motivated.  We  take: 

(86)  g(X)  =  aX 

In  this  case  ct  describes  the  strength  of  the  random  fluctuation  of  the  individual 
growth  rate  r(X). 

As  discussed  before,  (85)  is  equivalent  to  a  diffusion  equation 


(87) 
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where  p(x,  t)  is  the  probability  density  that  the  individual  number  X  shows  the 
value  X  at  time  t.  The  coefficients  are  given  by: 


(88)  B(x)  =  oV 

(89)  A(x)  =  xr(x)  =  xr(l-x/K) 


Next  we  show  that  the  population  dynamics  with  environmental  noise,  equation 
(87),  can  be  approximated  by  the  birth  and  death  process  (83).  We  first  discretize 
the  right-hand  side  of  (87)  by: 


(90) 

(91) 


d  j,/  ^  B{x  +  h)-B{x) 

—B{x)  = - ; - 

dx  h 


D/  \  B(.x  +  h)-2B{x)  + B{x-h) 
—Bix)  = - 


dx 


with  the  discretization  length  h  for  variable  x  (h  is  usually  set  to  1).  The  deriva¬ 
tives  of  A(x)  and  P(x)  are  handled  in  the  same  way.  A  simple  rearrangement  of 
(87)  with  (90)  and  (91)  and  the  definition: 

(92)  p(nh,  t)  =  p^(t) 


results  in  an  equation  (83)  of  the  birth  and  death  process.  The  birth  and  death 
rates  in  this  case  are: 


(93) 

II 

(94) 

11 

1  .Birth)  ^  A{nh) 
h  ■ 


h 


In  the  final  step  we  combine  both  discrete  models,  i.e.,  the  discrete  description  of 
environmental  noise  by  (83)  with  (93)  and  (94)  at  high  individual  ntimbers  and 
the  model  (83)  which  describes  demographic  noise  at  low  individual  numbers. 
Therefore,  our  final  stochastic  d3Tiamics  model,  which  is  valid  for  both  small  and 
large  population  communities,  is  the  birth-death  model  (83)  with  the  following 
birth  and  death  rates: 


(95) 


\  -  ^  [ctV  -I-  Rn  +  nr(n)] 
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(96)  M  +  Rn  -  nr(n)] 

”  2 

where  Rn  is  the  demographic  noise,  and  r(n)  is  the  density  dependent  intrinsic 
growth  rate.  The  simplest  form  of  r(n)  is  the  type  of  logistic  growth  rate,  i.e., 
r(n)=r(l-n/K). 

A  natural  community  usually  has  more  than  one  species.  The  Lotka-Volterra  or 
the  Kolmogorov  model  can  model  competition  among  species  if  there  are  only  two 
or  three  species  in  the  commimity.  Using  the  Lotka-Volterra  model  or  the  Kol¬ 
mogorov  model  in  a  highly  diverse  community  may  cause  chaos  because  of  the 
feedback.  This  problem  can  be  solved  by  introducing  the  growth  rate  function 
r(n)  as  some  other  form  of  the  resources.  For  example,  if  the  total  available  re¬ 
source  of  the  community  is  K,  and  the  used  resotmce  at  time  t,  is  K^,  then  we  can 
define  the  growth  rate  of  the  i-th  species  as  r'(n)=r'*(l-K/K)  or  some  other  form 
which  can  express  the  competition  relationship  among  species. 

For  the  combined  model,  we  want  to  emphasize  the  following  points.  If  the  den¬ 
sity  dependence  term  n/K  and  the  noise  term  a^n*  in  equations  (95)  and  (96)  are 
omitted,  we  get  the  model  with  only  demographic  noise.  This  model  is  motivated 
by  simple  biological  arguments  and  has  been  used  with  small  popiilations  by 
many  authors.  When  the  demographic  noise  term,  Rn,  is  omitted,  the  resulting 
model  is  a  diffusion  equation,  which  has  been  considered  by  many  authors  to  be 
a  description  of  environmental  noise. 
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4  Stochastic  Estimation  and  identification 


Quantitative  population  study  has  been  very  popular  for  several  decades  and  its 
models  have  been  well  developed.  These  models  fall  in  the  categories:  dynamic 
models  and  stochastic  models.  In  their  book,  Modeling  Fluctuating  Populations, 
Nisbet  and  Girniey  (1982)  present  the  most  common  types  of  dynamics  and  sto¬ 
chastic  models  and  the  application  conditions  for  each. 

One  of  the  basic  discrete  dynamics  models  of  plant  popcdation  is  the  logistic 
model  (Nisbet  and  Gurney,  1982): 

(1)  k  =  0, 1, ... 

where,  and  are  the  population  size  at  time  k-i-1  and  k,  respectively,  h(x,t)  is 
the  harvesting,  I  is  immigration  or  fertilization,  and  a  and  b  are  model  parame¬ 
ters. 

Consider  a  more  general  deterministic  system: 

(2)  =fk(Xk,uJ,  k  =  0, 1, ... 

where  ei?”  is  the  state  and  sR”  is  the  input  at  time  k. 

One  property  of  equation  (2)  is  that  if  the  current  state  x,^  and  the  input  sequence 
of  the  system,  u^,  u^^j, ...,  u^^„  are  given,  the  futxire  states,  x^^j,  x^.^2>  •••> 
determined  exactly.  That  is,  the  deterministic  model  (2)  makes  exact  prediction 
of  the  future  population  of  the  plant  system. 

Under  special  laboratory  conditions  and  in  some  isolated  environments,  popula¬ 
tion  growth  may  act  in  a  deterministic  way.  However,  most  natural  populations 
do  not  behave  as  nicely  as  described  by  the  deterministic  model  (Hallam  1986). 
Populations  fluctuate  about  some  deterministic  trend  due  to  demographic  sto- 
chasticity  and  environmental  stochasticity  (Turelh  1986).  In  Chapter  3  we  dis¬ 
cussed  how  to  model  stochastic  systems. 

A  stochastic  model  is  able  to  predict  the  probability  that  at  a  given  time,  the 
population  will  be  of  a  particular  size.  Given  the  probability  or  the  conditional 
probability  of  future  behavior,  the  expectation  of  the  population  size  can  be 
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determined  as  it  is  in  the  deterministic  model.  In  the  real  world,  however,  the 
probability  distributions  of  most  plant  populations  are  rarely  known.  Because 
we  often  do  not  have  descriptive  information  for  the  probability  distribution,  past 
data  must  be  examined.  The  availability  of  observations  from  the  past  allows  for 
an  alternative  method  of  population  modeling:  the  technique  of  regression.  This 
technique  has  become  a  common  tool  in  estimation  and  prediction  in  many  areas. 

One  serious  problem  with  regression  models  is  the  assumption  of  constant  coeffi¬ 
cients.  Because  of  this  assimiption,  underlying  phenomena  can  never  be  de¬ 
tected.  Also,  the  reliability  of  the  estimate  of  regression  depends  on  the  sample 
size.  A  regression  model  usually  requires  a  minimal  sample  size  of  30.  However, 
population  data  that  spans  30  or  more  years  is  not  common.  Sometimes  we  can¬ 
not  afford  to  have  a  large  sample  size.  Moreover,  observed  data  used  in  regres¬ 
sion  is  assumed  to  be  correct.  This  is  not  always  the  case.  In  the  real  world, 
most  natural  communities  are  full  of  noise.  Observed  data  contains  errors  such 
as  system  noise  and  measurement  error.  For  some  regression  models,  these  er¬ 
rors  not  only  create  a  large  variance,  but  also  produce  a  biased  estimate  (Gert- 
ner,  Cao,  and  Zhu  1995). 

To  get  a  better  grip  on  these  problems,  we  introduce  the  theory  of  optimal  esti¬ 
mation.  An  estimator  is  a  process  by  which  information  can  be  extracted  from 
data,  i.e.,  to  infer  desired  information  by  filtering  out  the  noise  from  the  data. 
Because  the  estimator  combines  the  descriptive  information  amd  the  data  infor¬ 
mation  of  a  system  to  form  an  estimate,  the  estimate  usually  has  a  lower  vari¬ 
ance  than  a  conventional  regression  estimate.  The  most  common  estimator  is 
the  Kalman  filter  (Kalman  1960).  The  filter  describes  how  to  process  the  meas- 
mement  data  for  a  given  linear  system.  The  theory  of  optimal  application  has 
been  successfully  applied  in  a  broad  range  of  areas.  These  areas  include  signal 
processing  in  communications,  tracer  studies  in  medicine,  statistical  image  en¬ 
hancement,  estimation  of  traffic  densities,  chemical  process  control,  satellite  or¬ 
bit  estimation,  unclear  reactor  parameter  identification  (Gelb  1974),  and  dendro- 
cHmatology  (Visser  and  Molenaar  1988). 


Optimal  Estimation 

An  optimal  estimator  is  a  computational  algorithm  that  processes  measurements 
in  order  to  deduce  a  minimum  error  estimate  for  the  state  of  a  system.  It  does 
this  by  using  knowledge  of  system  and  measurement  dynamics,  assumed 
statistics  of  system  noise  and  measurement  errors,  and  with  initial  condition 
information.  The  advantages  of  this  type  of  data  processor  are  that  it  minimizes 
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the  estimation  error  in  a  well-defined  statistical  sense  and  that  it  uses  all 
measinement  data  plus  prior  knowledge  about  the  system. 

There  are  three  types  of  estimation  problems:  filtering,  smoothing,  and  predic¬ 
tion  (Gelb  1974).  When  the  time  at  which  an  estimate  is  desired  occiu's  at  the 
last  measurement  point,  the  problem  is  referred  to  as  filtering.  When  the  time  of 
interest  falls  within  the  span  of  available  measurement  data,  the  problem  is 
termed  smoothing.  When  the  time  of  interest  occurs  after  the  last  available 
measurement,  the  problem  is  called  prediction.  In  the  following  sections  we  will 
discuss  linear  filtering  and  prediction,  and  two  techniques  of  nonlinear  approxi¬ 
mation.  For  convenience,  we  consider  a  stochastic  system  without  input.  Incor¬ 
porating  a  known  input  is  straightforward. 

Optimal  Linear  Filtering 

First  let  us  consider  a  discrete  stochastic  system  without  input  control: 


(3) 


+  Gk^k 

yk  =  c,,x^+^kVk’ 


where  eR\y^  eR'’,w^  K’  are  possible  time- 

varying,  known  matrices  of  appropriate  dimension,  x  and  y  are,  respectively,  the 
state  space  and  observation  space.  The  basic  random  variables  {x^,,  w^,,  ...,  V(,,  ...} 
are  all  independent  and  Gaussian,  with  Xq  ~  iV^(0,S(,),  w^~N(0,Q), 

~  N(0,R).  The  covariances  are  all  known.  The  available  information  at  time 
k  is  z*  =y:=(yi,y*_i,—,yo)-  The  random  variable  x^,  x^^j,  and  y'‘  are  jointly 
Gaussian;  denote 


PkikMy')- 

Pk+\\k  ('^i+1  \y  )  ~ 


By  definition, 

(5)  , 

'^k\k■=E{ix^-x^^^)(x,-x^^,^)  \y  }. 


Similarly, 


(6) 


To  obtain  the  recursion  rule,  we  take  the  following  steps. 


Step  1 

From  (3), 

(7) 

For  convenience  denote 

^k+]\k  •  ~  ^k+i  ~  » 

(8) 

^k\k'~  ~^k\k- 

From  (3)  we  have 

(9)  •^i+||i-  =  ^k^klk  "^^k^k 

By  the  independence  of  and  w^, 

(10)  2:,„,  =  4i,,^;+G,eG; 

Step  2 

Denote 

yk\k-r=E{y;^\y'‘-'},and 

yk\k-V~yk  ~yk\k-\ 


Since  v^  and  y‘‘‘^  are  independent,  this  gives 

(12) 

yk\k-\  ~  ^k^k\k-\  "^^k^k- 

By  the  independence  of  and  v,^ 


(13)  1^,.,:=  cov(y,,_,)  =  + 

From  (12)  and  the  independence  of  and 

(14)  ^^kyk\k-\  ~  '^k\k-\^k 


T 


,  we  also 


42 


CERL  TR  00-5 


Step  3 

By  the  innovation  principle,  we  have 
(15)  ^k\k~^{^k\y 


=  £  {;c*  I  y*-’ }  +  {Ex,  )"  y*,* 


“  f^k^^k  ]  7kVi-\- 


and 


(16) 


Therefore,  the  conditional  density  can  be  obtained  from  the 

following  recursion  relations. 

(17) 


^0|0 


^A+1|A+1  A'+lQ+l)^ 


A+ljA, 


^A+lIA  ~  ^k^k\k^k  '^^kQ^k 


2^010““  (7-IoCo)Z  Oj 


here 


A*  —  +  ^k^^k  ] 

4=2„C[q,2oC+^4R«»’']-' 

The  recTxrsion  algorithm  (17)  specifies  the  transition  function  of  the  information 
state.  It  is  known  as  the  discrete  Kalman  filter.  The  matrix  L  in  the  recursion  is 
called  the  Kalman  gain  matrix  (Gelb  1974).  It  can  be  shown  (e.g.,  Otter  1978) 
that  the  ordinary  least  squares  (OLS)  fitting  procedure  is  a  special  case  of  the 
Kalman  filter  (17). 

The  transition  from  the  discrete  to  the  continuous  formulation  of  the  Kalman 
filter  is  straightforward  by  setting  time  k  and  k-^l  in  the  discrete  case  to  time  t^ 
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aBd  respectively,  and  letting  =  A/  ->  0.  First,  consider  the  following 

continues  stochastic  system: 

(18)  x  =  Fx  +  Gw 

y  =  Cx  +  v 


where  w,  v  are  white  noise  processes  with  mean  0  and  spectral  density  matrix  Q 
and  R,  respectively.  In  this  continuous  case,  v  is  equivailent  to  H^v^  in  the  dis¬ 
crete  case.  To  make  the  transition,  let  us  first  rewrite  the  discrete  system  as  a 
difference  equation. 

(19)  X,  =Aj^x, 


± 


(A-O  y 
A/  ■*/,. 


From  this  difference  we  have  the  following  equivalence,  valid  in  the  hmit  as 
=  Ai->0: 

A,^I  +  FM 

In  the  discrete  case  we  have  shown  that: 

(20)  =  + 


This  is  now  rewritten  as: 

Af^  — ^  I  "i-  FAi 

In  the  discrete  case  we  have  shown  that: 

(21)  =  (/-fFAi)Z^n.(/-FFAi)[  -f  G/^QG^At 


Expansion  3rields: 

(22)  ^*+11*  ^*1* 


+  (FE,,  -H  E„,F"  +  G,QGl  )At  +  0(A?^) 


From  the  discrete  model  we  also  have: 


(23) 


Insert  (23)  into  (22)  and  rearrange  the  terms  we  have: 


(24) 

-F4QI.^_,  -4QI.„_,F’'  +  0(A/) 

Now  let  us  exam  the  term  /  A/ 

'h^-'k  -'b^k\k-\^k  [Q^*|Ar-lQ  ^k^k  ] 

(25)  =  \.^k^k\k-\^k  ^  ■*“  ^k^^k 

=  2:.»-,Q''[Q2.MC>+flr' 

Taking  the  limit,  we  get: 

(26)  limiZ,  =SC^i?'' 

Ar->0 


and, 

(27) 


t{t)  =  F(02(0  +  nt)F^it)  +  G{t)Q{t)G^{t) 

-l.{t)C^{t)R-\t)Cit)m 


2(0)  =  So 

By  similar  manipulation,  we  have  the  state  estimate  of  the  continues  system. 

(28)  kt)  =  +  L{t)[y{t)  -  C(0x(t)] 

The  matrix  gain  is: 

\  nt)C^(t)R-'(t),  £[w(/)v'(t)]  =  0 

(29)  «''-|[£(,)c'-(,)  +  C(<)4.(()]F-W,  £(M-(0v'(T)]  =  *wa(/-T) 


Optimal  Prediction 

The  prediction  model  can  be  derived  from  the  filtering  model  in  a  straightfor¬ 
ward  fashion.  The  one-step  predictor  is  given  by  equations  (7)  and  (10). 
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Nonlinear  Estimation 

Now  let  us  consider  a  more  general  case  described  by  the  nonlinear  Stochastic 
Difference  Equation  (SDE)  with  discrete  observations. 

(33)  X  =  f(x,t)  + W(t) 

yk  k=i,2, ... 

where  f(x,t)  and  h(x,t)  are  nonlinear  functions,  w(t)  is  Gaussian  noise  with  mean 
0  and  having  spectral  density  matrix  Q(t),  and  ~  N(0,  ) .  One  further  com¬ 

plicated  mode  is  the  form  of  x  =  f(x,t)  +  g(x,t)w(t).  In  this  case  a  theory  for  es¬ 
timating  x(t)  cannot  be  developed  within  the  traditional  framework  of  mean 
square  stochastic  calculus  because  the  r.h.s.  of  the  equation  is  not  integratable  in 
the  mean  square  sense.  This  difficulty  is  overcome  by  formulating  the  nonlinear 
filtering  problem  within  the  context  of  Ito  calcffius. 

We  first  discuss  the  nonlinear  system  with  the  form  of  (33).  There  are  two 
widely  used  hnearization  techniques:  tnmcated  Taylor  expansion  and  statistical 
approximation  (Gelb  1974). 

The  Taylor  expansion  method  is  to  write  nonlinear  function  f(x,t)  and  h(x,t)  as  a 
Taylor  expansion  about  the  current  estimate  of  the  state  vector. 
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f(x,t)  =  f(x,t)+fU(x-x)+... 
h(x,t)  =  h(>(,t)  +  f|,^,(x-x)+... 


Discarding  all  the  terms  with  second  or  higher  order,  we  obtain  a  linear  ap¬ 
proximation  of  f(x,t)  and  h(x,t)  about  the  current  estimate  of  the  state  vector. 
The  derivation  of  the  linear  fQter  has  been  discussed  in  the  previous  section. 

To  obtain  an  accurate  estimate,  we  need  higher-order  filters.  One  method  by 
which  the  estimate  can  be  improved  is  by  repeatedly  calculating  ,  Lj^,  and 
each  time  hnearizating  about  the  most  recent  estimate.  Another  method  is 

to  include  more  terms  in  the  expansions  for  f(x,t)  and  h(x,t). 

The  statistical  linearization  is  to  seek  a  linear  approximation  for  a  vector  func¬ 
tion  fix)  of  a  vector  random  variable  x,  having  probability  density  function  p(x). 
The  approximation  is  made  by  change  functions  fix)  and  h(x)  in  an  approximate 
linear  form. 

/(x)s£[/(x)]+iV/x-x) 

(35/ 

h,(x)^E[h,(x)]  +  N,(x-x) 

where  Nj  and  N,^  are  called  the  function  gain  matrices.  They  are  estimated  by 
using  the  technique  of  minimtun  mean  square  error.  These  estimates  are  given 
by: 


(36)  ^  ^  ^ 

Where  E(fa^),  f,  x  are  expectations  calculated  assvuned  x  ~  iV(x,S). 

are  expectations  calculated  assumed 

x(^)  ~ 

Because  nonlinear  functions  fix)  and  h(x)  are  linearized,  the  approximate  opti¬ 
mal  filter  for  the  linearized  system  now  can  be  found  as  discussed  above.  In 
many  instances,  the  statistical  linearization  technique  has  better  performance 
than  the  Taylor  expansion  method.  However,  the  decision  as  to  which  types  of 
filters  should  be  used  in  a  particular  application  depends  upon  their  computa¬ 
tional  complexity  and  relative  performance  as  observed  from  realistic  computer 
simulations. 
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Let  US  consider  a  more  general  nonlinear  stochastic  system. 

(37)  x,=f{x„t)  +  G{x„t)w^ 

yh 

Now  suppose  that  we  generate  a  reference  deterministic  trajectory  xit),  with 
given  x(io),  satisfying: 

(38)  j  =  /(J(0,0 
Define: 

(39)  &Cj=x, -3c(i), 

=  y,  -  yW,  and 

Linearizing  (37)  about  x(t)  by  a  Taylor  series  expansion,  we  obtained  the  line¬ 
arized  discrete  system: 

(40)  ^ 

8yi^=M[t^,x(tk)]^t,+Vk 

Where  3c (i;^)]  is  the  state  transition  matrix.  That  is, 

(41)  doit,  x)l  dt  =  f  it)Oit,  t)  ,  and 
0(^,x)0(x,^)  =  a)(^,^) 

The  new  noise  term,  ^ ,  is  given  by: 

(42)  _  =J  a)(/^^,,T)G(T)Jw^ 


which  is  a  form  of  Ito’s  stochastic  integral.  From  the  theorem  of  Ito’s  stochastic 
integral,  it  is  easy  to  see  that  }  is  a  white  Gaussian  sequence  with  a  mean  of 

0  (Jazwinski  1970),  {w^^  }  ~  N\0,Qik  -f  1)],  where 
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(43)  e(t + 1)  =  r  *(<..,  i)G(i)e(T)G"(t)«>’‘((,,,,T)<ir 

Applying  the  linear  filter  to  the  linearized  system  (40),  we  obtain  the  extended 
Kalman  filter  for  the  nonhnear  system  (37)  (Jazwinski  1970). 

~  I,.  ■ 

~  ^^h+v  3 

The  Kalman  gain  is: 

(46)  [4+i> ^ 

To  improve  the  reference  trajectory,  we  need  some  iteration  algorithms  in  which 
the  estimate  Ti,.,!  =  “ ^il} 

can  be  improved  by  repeatedly  calculating  i=  1,...,/,  L]j,  and  each  time 

linearizating  about  the  most  recent  estimate.  The  following  is  the  iterated  ex¬ 
tended  Kalman  filter  with  (45)  replaced  by  the  iterator  (Jazwinski  1970): 
nM  =  --ni]}.  with 

i  =  l,...,l  and  =  i)/-  The  iteration  starts  with  T|j  =  ,  and  terminates 

when  there  is  no  significant  difference  between  consecutive  iterates. 


(44) 


(45) 


Stochastic  Identification 

A  Kalman  filter  can  yield  optimal  performance  for  a  linear  stochastic  system. 
Consider  the  discrete  stochastic  system  (3) 

(47)  +  B^Uk  + 

yk=^kH+tik^k, 
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The  Kalman  estimator  (conditional  mean  tmder  given  information 

is  given  by: 

(48)  ^A+uft+i  “  ^k+^yk+\~^k+i^k^k\k^’ 

The  basic  assiomption  of  the  Kalman  estimator  is  that  the  descriptions  of  A,  B,  G, 
C,  H,  Q,  R,  and  Eq  [as  defined  in  (3)]  are  correct.  As  a  practical  fact,  this  is  usu¬ 
ally  impossible. 

Estimations  of  these  quantities  must  be  made  and  improved  by  the  information 
fi”om  observation.  Therefore,  we  need  to  develop  some  procedure  that  provides 
the  best  descriptions  of  A,  G,  C,  H,  Q,  R,  and  Eq  from  observations  and  other 
prior  knowledge.  This  problem  is  called  parameter  identification. 

Consider  the  following  system: 

p  h  r 

(49)  y,  =  ^  -l-  ^  c, w,_, 

1=1  /=1  1=0 

This  is  an  Autoregressive  Moving  Average  Model  (ARMAX)  model.  The  first 
term  on  the  r.h.s.  of  equation  (49)  expressing  the  dependence  of  the  current  out¬ 
put  on  its  own  past  values  is  the  autoregressive,  or  AR,  term.  The  second  term  is 
the  external  inputs  (or  control  inputs)  of  the  system.  The  last  term,  which  is 
called  the  moving  average  (MA),  is  a  moving  combination  of  independent  random 
variables  v . 

k 

Using  a  shift  operator,  the  ARMAX  model  (49)  can  be  expressed  as: 

(50)  A(q'')yj^  =  +  C(q~')w^ 

where  q~^  is  a  backward  shift  operator,  i.e.,  q~'y^  =  y^^i,  and 

A{q~')  =  \-¥a^q~'  +--^a„q~’’ 

(51)  B(^q  ^)  =  bQ+b^q  ^ +-“  +  b„q  * 

The  stationary  condition  of  an  autoregressive  process,  such  as 

-I- — is  that  all  roots  z,  of  its  pol5momial,  1-<3|Z - 

have  modioles  |zo|>  1.  Also,  a  moving  average  process  of  order  r,  MA(r),  can  be 

written  as  an  infinite  AR  process  if  all  the  roots  of  its  polynomial, 
1  -  c,z - c^z'^ ,  have  modules  greater  than  1.  Such  an  MA  process  is  called  in- 
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vertible.  Without  special  notice  we  henceforth  assume  that  the  MA  term  in  the 
AEMAX  model  is  invertible. 

Least  Square  Mean  Estimates  (LSE) 

Suppose  we  have  the  data  and  we  beheve  that  the  follow¬ 

ing  model  fits  the  data. 

(52) 

where  Wy.  is  the  error  term  with  a  mean  of  0. 

Let  F  be  a  s-field  generated  by  the  data  ,  •  •  •,  ,  y, ,  •  •  •,  } .  Then,  {w,^,  k  £  n} 

n 

is  F^-measurable  because  can  be  deduced  from  yj^  and  Fj^  through  the  relation 

We  also  assume  that  {w^}  is  a  Martingale  difference  sequence  with  respect  to  the 
increasing  sequence  of  s-field  {3  J  which  satisfies 

E[wl\^^_^]  =  (7^,  a.s.,  for  all  k 

and 

SUp.£’[|Wt|“|J,_,]<oo,  a.s.,  for  some  a  >  2. 

k 

The  least  square  method  is  to  choose  0  to  minimize 

(53) 

k=0 

By  setting  =  0 ,  we  obtain  the  LSE  of  6 

(54) 

*=0  t=0 

Suppose  one  more  datum  becomes  available;  then  we  can  obtain  the 

new  LSE  from  the  old  LSE0„ . 


(55) 
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n 

Where  •  To  verify  (55),  we  substitute  the  equation 

*=o 


k  =  by  =  K-xK  and  write 

Jt=0  /b=0  ^=0 


-o  o 

n 


T 

n 


Let  P  =  i? We  obtain  the  LSE  recursions. 


(56) 


P  =P  , 

n  n—\ 


l40^p  o  ^n-1  ^n^nr 


'n'^n-X'^n 


0*^  ^  ^  " 

For  convenience,  we  define  9^  =  9^  -  9^,  and  -  0^,0^ 


Theorem  (LSE  convergence) 

Consider  the  LSE  algorithm  defined  by  (55)  eind  applied  to  model  (52),  then 


(i)  £(!»,-»,)' =0(logrr(V,))  a.s. 


k=0 


(ii)  = 


n-] 

(ni)  ^  ^  bare  i=l, ...,  2p+2, 


n->oo 


Moreover,  if  =  +°o  and  >  e/  for  all  large  n,  then 


lim^n  =  ^  a.s. 

n-+oo 

Proof 

We  will  not  prove  results  (i)  to  (ii)  here.  They  can  be  proved  by  using  Martingale 
convergence  theorem  and  Kronecker’s  lemma  (i.e.,  two  real  valued  sequences  {xj 

00  N 

and  {r  }  satisfying  r  >  0,  linir*  =  ^  ^^en  liniTv 

*->®  *=1  *  N-t<B  "  i=l 

give  the  proof  of  (iv),  the  convergence  of  the  estimate. 
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Since 


n~\ 

1  +  ^  (O'^.  )^  <  1  +  Tr ^  ,  from  restilt  (iii)  we  have 

*=0  i=0 


lim 

/7->CO 


1 

n-l  - 

Trl<PtOl 

*=0 


*=o 


and 

k=0 


n-l 


7*1-1 


Tr  t=0 

A=0 


x[ 


n-\ 


TrZ<t>i0[  k=0 
i=0 


Taking  the  limit  of  the  equation  above,  we  obtain  result  (iv). 


Suppose  we  have  the  data  one-step  predic¬ 

tion  made  by  0„  is 


Denote  0;  is  the  i-th  element  of  0,  then  the  feedback  law  is  explicitly  given  by 

“n  =  i::S^A(^)yn +-+K 

(57)  +^„(^ +  2)w„_, 

-yn*^-K(Q+h  +  l)y„ - K(9  +  f^  +  '•)  A-r+I 

Similar  results  can  be  obtained  for  stochastic  gradient  (SG)  algorithm.  The  SG 
algorithm  is  regarded  as  a  simplification  of  the  LSE  recursions  since  the  scalar 
gain  in  the  SG  algorithm  is  only  the  trace  of  the  matrix  gain  in  the  LSE  algo¬ 
rithm. 

Conditional  Mean  Estimates 

Let  us  rewrite  the  ARMAX  model  (49)  as 
yk+\  ~  ‘I’*  ^  ^*+1 

where  =  (y,,---,y, =  . 

The  true  parameter  6°  is  xmknown  and  we  want  to  identify  it. 
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Theorem  (conditional  mean  estimate) 

Assume  that  6^  ~  N{G,  E),  {y^}  is  determined  by  =®[0°  +  W;t+i,  {W;}  is 
identical  independent  distribution  (iid)  with  ~  N{0,O^),  and  {y^,^^,s<k}  is 
independent  of  {w^,s>k  +  l} .  Then  ,y^^“‘,y„}  minimizes 

and  0„  and 

P„:=  E{(eP-  G„)(eP  -  satisfy  the  recursions 


(58) 


P„  =  - 


This  theorem  says  if  {w^, }  is  Gaussian  white  noise,  the  recursions  for  the  condi¬ 
tional  mean  of  0®  xmder  Bayesian  formulation  coincides  with  the  LSE. 


Proof 

Let  F(e)=E{||0°  -  9if  |Oo,-,(D„_,  ,y„-,y„} 

minr(e)}=^^=« 

e 

^2E{\ef-ei  •,y.)=0 

=>e=£{e»l<i>„,.-,(i>,^, y,} 

To  prove  the  recursion  (58),  we  rewrite  the  ARMAX  model  as  a  state  space  repre¬ 
sentation. 

G,-N(G,  E), 

>'*+1  + 

Then  we  can  use  the  Kalman  filter  as  discussed  in  the  previous  section 

^k+\\k-i-\  “  ^k^k\k  ^k+]\.yk+]  ^  ^k-k\^k^k\k\ 

^A:+1|A-+1  ~  '^A'+lQ+l)^A-+ltA', 
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and  make  the  following  exchanges 

^k\k  ^  ^k 
^  ^k 

A/^oI,  Gj  <te>0,Q  <»0[, 


We  obtain  the  recursion  (58). 
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5  STELLA  II  Modeling  Process 


In  the  previovis  chapters,  we  have  chosen  and  defined  theoretical  statistical 
measxires  of  plant  diversity  and  developed  theoretical  models  for  plant  communi¬ 
ties.  The  stochastic  dynamics  models  developed  contain  two  main  components: 
deterministic  process  and  stochastic  process.  The  deterministic  model  controls 
the  biological  dynamics  of  the  plant  commvmities,  while  the  stochastic  process 
simulates  the  biological  and  ecological  fluctuations.  When  we  have  a  good  un¬ 
derstanding  of  the  dynamics  and  stochastic  behaviors  of  the  plant  commionity  we 
are  studying,  the  stochastic  dynamics  model  is  able  to  make  respectable  predic¬ 
tions  for  the  community.  Once  we  have  chosen  to  use  a  stochastic  model,  there 
remains  the  task  of  selecting  the  mathematical  approach  to  follow  in  its  analysis. 
Figure  1  illustrates  the  sequence  of  decisions  involved  in  choosing  an  appropriate 
model.  The  first  decision  is  the  specification  of  birth  and  death  probabilities. 
The  second  is  to  determine  the  size  of  fluctuation  due  to  demographic  variation. 
Finally,  a  choice  of  representation  of  environmental  variation  must  be  made. 
Following  the  flow  of  the  diagram,  we  can  select  one  of  the  mathematical  models 
discxissed  in  Chapter  3  as  om  simulation  model. 


Figure  1.  Sequence  of  choosing  a  stochastic  model. 


As  we  pointed  out  before,  a  stochastic  model  is  able  to  predict  the  probability 
that  the  population  wiU  be  of  a  particular  size.  Given  the  probability  or  the 
conditional  probability  of  future  behavior,  an  estimate  of  the  poprdation  can  be 
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determined  as  in  the  deterministic  model.  In  the  real  world,  however,  we  usually 
do  not  know  the  probability  distribution,  and  we  do  not  know  the  details  of 
population  dynamics.  Instead,  we  must  look  at  past  data.  In  this  case,  estima¬ 
tion  models  discussed  in  Chapter  4  will  be  our  choice  of  simulation  models. 

U.S.  military  installations  across  the  country  cover  thousands  of  plant  commu¬ 
nity  types.  With  the  development  of  the  U.S.  Army  Land  Condition  Trend 
Analysis  (LCTA)  program.  Army-wide  inventory  plots  have  been  established  and 
field  data  have  been  collected  on  many  Army  installations  since  1989.  For  this 
reason,  we  choose  the  estimation  model  presented  in  Chapter  4  as  our  major 
STELLA  simulation  model.  (STELLA  is  a  software  modeling  tool.)  The  dynam¬ 
ics  model  is  used  only  for  testing  purposes. 

Consider  the  discrete  stochastic  system  presented  in  Chapter  4. 

(1) 

yk  ~  ^k^k  ’ 

where  Xf^  GR'^,yi^  &R^,Vi^  ei?*;  A^,  G^,  C^,  and  are  possible  time- 

var3dng,  known  matrices  of  appropriate  dimension,  x  and  y  are,  respectively,  the 
state  space  and  observation  space.  The  basic  random  variables  {x^,  w^,,  ..,  v^,  ...} 
are  all  independent  and  Gaussian,  with  Xq  ^  N ,  w^~'N(0,Q), 

Vi^  ~  N  (0,  R) .  The  covariances  are  all  known.  The  available  information  at  time 
k  is  The  random  variable  x^,  x^^j,  and  y""  are  jointly 

Gaussian.  As  derived  in  Chapter  4,  the  recmsion  scheme  was  obtained  by  Kal¬ 
man  filtering: 

(2)  =  ^k^klk 


^010  ~  ^oyo, 

^A+Ui+1  ~  “  ^k*fik+l^^k+-ilk, 

So,o  =  (/-4Co)So, 


1 


Here 
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This  recursion  specifies  the  transition  fiinction  of  the  information  state.  It  is 
known  as  the  discrete  Kalman  filter.  The  matrix  L  in  the  recursion  is  called  the 
Kalman  gain  matrix.  The  following  STELLA  model  uses  this  recursion  scheme 
to  model  the  changes  of  the  biodiversity  and  population  of  plant  communities. 


The  STELLA  Model 

In  this  study,  the  STELLA  model  is  built  based  on  the  recursion  scheme  (2). 
Figure  2  is  a  STELLA  map,  which  shows  the  fo\ir  basic  parts  of  our  STELLA 
model:  global  variables,  local  variables,  the  Kalman  filter,  and  the  outputs.  A 
variable  is  called  global  if  it  is  determined  by  the  community  type  and  is  inde¬ 
pendent  of  the  species  type,  such  as  the  Shannon  index,  total  number  of  species, 
total  population,  total  number  of  observations,  error  hmit,  natural  input,  and 
human  input  (these  are  discussed  in  the  following  section).  A  local  variable  var¬ 
ies  with  the  type  of  species.  Local  variables  include  species  abundance,  species 
tolerance  to  environmental  changes,  and  model  parameters.  The  model  of  the 
Kalman  filter  is  based  on  species.  Although  the  structure  of  the  model  is  the 
same  for  all  species,  the  inputs  and  outputs  are  different  for  different  species. 
The  outputs  include  all  the  major  results  of  the  simulation  both  in  table  and 
graphic  format. 

Figure  3  is  a  STELLA  diagram  that  shows  the  structure  of  the  estimation  model 
of  the  Kalman  filter.  The  parameters  emd  components  in  the  model  are  defined 
in  the  recursion  scheme  (2).  This  model  simulates  the  changes  of  the  population 
of  a  single  species.  A  multispecies  community  needs  multiple  copies  of  this 
model  in  which  the  structure  is  the  same  but  the  parameters  are  set  differently. 


Figure  2.  The  STELLA  map  of  flow  structure. 
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STATE 

NEW  STATE  i - 1  OLD  STATE 


VARIANCE 

NEW  VAR  I - 1  OLD  VAR 


CURR^XOBS 


KALMAN  GAIN 


OBSERVATION 


KALMAN  GAIN 


Figure  3.  The  STELLA  structure  of  the  Kalman  filter. 
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Setting  of  the  STELLA  model 
Global  Variables 

Number  of  observations:  is  the  number  of  repeated  times  of  observation. 

Error  limit:  is  the  tolerance  of  inventory  errors  that  include  sampling  and  non¬ 
sampling  errors. 

Human  input:  is  the  quantified  hiaman  activities,  including  the  destruction  and 
improvement  activities  done  to  the  community. 

Natural  input:  is  scaled  environmental  variation  and  natural  catastrophe.  The 
growth  of  a  plant  population  is  determined  not  only  by  its  intrinsic  reproduction 
and  survival  capabilities  but  also  by  its  surrounding  climatic,  topographic,  and 
geologic  conditions.  Some  of  these  environmental  factors  have  httle  variation 
over  time,  such  as  topography  and  soil,  while  others,  such  as  temperature  and 
precipitation,  change  daily.  Some  of  these  changes  are  important  to  the  popula¬ 
tion  growth,  and  some  of  them  are  not.  We  choose  those  factors  that  are  impor¬ 
tant  to  the  community  and  change  over  a  relatively  short  period  of  time.  Catas¬ 
trophes  can  be  thought  of  as  the  extremes  of  environmental  variation.  These  are 
events  that  affect  either  reproduction  or  survival.  Catastrophes  include  habitat 
destruction,  flood,  fire,  disease,  drought,  storm,  etc.  We  may  be  able  to  define  the 
impacts  of  these  catastrophes  on  the  community  of  interest  by  examining  the 
historical  records  of  catastrophes  and  popiilation  changes  of  the  community. 

Total  population:  is  the  sum  of  the  populations  of  all  species  in  the  community. 

Total  number  of  species:  is  the  total  niimber  of  species  appearing  in  the  commu¬ 
nity  during  the  simulation. 

s 

Shannon  index:  is  defined  as  -^Pi  log(p,  ),  where  S  is  the  total  number  of  spe- 

/=1 

cies  and  p,-  is  relative  population  of  the  i-th  species.  The  Shannon  index  is  de¬ 
termined  by  both  the  richness  (number  of  species  of  the  commtmity)  and  the 
evenness  (relative  abimdance  of  each  species).  Given  the  total  number  of  spe¬ 
cies,  the  Shannon  index  reaches  its  maximxim  when  all  species  have  the  same 
abimdance.  That  is,  Max  (H)  =  log(5).  Comparing  the  Max  (H)  and  the  true 

value  of  H,  we  know  the  level  of  the  diversity  of  the  community. 

We  also  calculate  the  population  sizes  of  the  most  abundant  species  and  their 
variances.  In  most  natured  communities,  only  a  few  species  contribute  80 
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percent  or  more  to  the  total  population.  These  abundant  species  usually 
determine  the  quality  and  function  the  community  reaches. 

Local  Variables 

Model  parameters:  model  parameters  A,  G,  Q,  C,  H,  and  R  are  defined  in  the 
schema  (20).  In  our  study,  parameters  G,  C,  and  H  are  set  to  be  1.  A  is  the  net 
growth  rate.  When  a  community  is  in  a  stable  stage,  parameter  A  should  be  1. 
Thus  A  is  initialized  to  be  1.  Q  is  the  system  noise  term.  It  is  usually  unknown. 
As  we  proved  in  Chapter  4,  Q  does  not  affect  the  trend  of  estimates.  Even  when 
Q  is  initialized  large,  it  drops  very  fast  as  more  and  more  data  enter  into  the 
model.  R  is  the  combined  inventory  error,  which  is  determined  by  the  population 
sized  and  the  error  limit. 

Species  tolerance:  the  tolerance  of  a  species  is  the  capability  of  the  species  to  re¬ 
sist  environmental  changes.  It  is  one  of  the  intrinsic  characters  of  the  species. 
Different  species  have  different  sensitivities  to  changes  in  their  environment. 
Thus,  the  tolerance  index  is  determined  by  our  imderstanding  of  the  species. 

Observation:  the  observation  is  a  set  of  time-series  data  collected  from  the  field. 
The  size  of  the  data  set  affects  the  quality  of  the  simulation  results.  The  more 
data  we  have,  the  more  actuate  the  results  we  obtain  from  the  simulation. 

Kalman  filter 

Stocks,  Flows,  Converters,  and  Connectors  are  the  generic  building  blocks  of  the 
STELLA  language.  Stocks,  flows,  or  converters  represent  variables  and  parame¬ 
ters.  Connectors  link  related  variables  or  parameters.  This  section  provides  a 
brief  description  of  these  building  blocks,  and  shows  how  the  blocks  are  used  to 
build  the  Kalman  filter  (2).  STELLA  manual  (STELLA  II  1994)  provides  details 
of  modeling  in  STELLA. 

STOCK 

!•  In  our  model,  stocks  act  as  buffers  within  the  system  of  plant  popula¬ 
tions.  They  build  up  or  decline  whenever  their  associated  rates  of  inflow  and 
outflow  are  out  of  balance  with  one  another.  This  buffering  property  of  stocks 
leads  to  dynamic  simulation  of  the  system.  There  are  two  stock  variables  in  the 
Kalman  filter:  population  and  its  variance.  The  population  is  the  optimed  esti¬ 
mate  of  the  true  poptdation  of  a  given  species.  The  variance  is  the  minimized 
variance  of  the  population.  These  two  stock  variables  are  calcvilated  by  the  Kal¬ 
man  filter  (2). 


W  :  In  STELLA,  flows  consist  of  a  pipe  (or  conduit),  flow  regxdator,  and  ar¬ 
row.  Things  flow  through  the  conduit  in  the  direction  indicated  by  the  arrow. 
The  specific  volume  of  the  flow  is  calculated  by  the  algebraic  expression  in  the 
regulator.  In  the  Kalman  filter,  there  are  two  types  of  flows:  flow-in  and  flow- 
out.  A  flow-in  and  flow-out  connect  each  of  the  stock  variables  in  our  model. 
When  a  new  observation  comes  into  the  model,  the  regulator  in  the  flow-in  cal¬ 
culates  the  new  population  size  based  on  the  Kalman  filter,  and  adds  the  result 
into  the  stock.  Meanwhile,  the  flow-out  removes  the  old  population  from  the 
stock.  Therefore,  the  population  size  stored  in  the  stock  is  updated  as  the  recm- 
sion  scheme  (2)  shows. 

CONVERTER 

O  :  Circles  in  STELLA  represent  converters;  they  are  the  containers  for  all 
types  of  information  or  material  quantities.  As  their  name  implies,  converters 
transform  inputs  into  outputs  based  on  the  expressions  in  the  circles.  Unhke 
stocks,  converters  do  not  accumulate  flows  and  have  no  memory.  In  our  model, 
converters  represent  all  of  the  global  and  local  variables. 

Outputs 

The  outputs  of  the  model  include  total  population  size  of  the  community,  total 
number  of  species,  the  Shannon  index,  and  the  populations  of  the  five  most 
abimdant  species  of  the  commxmity  and  their  standard  deviations.  These  out¬ 
puts  are  presented  in  table  and  graphic  formats. 
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6  A  Case  Study:  White  Sands  Installation 

White  Sands  Missile  Range  was  selected  as  the  case  study.  The  White  Sands  in¬ 
stallation  is  located  in  the  desert  in  New  Mexico.  Although  the  species  composi¬ 
tion  here  is  not  complicated,  White  Sands  has  a  wide  range  of  land  covers,  from 
pru-e  sand  with  little  or  no  vegetative  cover,  to  highly  dense  plant  communities. 
A  single  model  may  not  be  able  to  model  the  whole  installation.  Therefore,  we 
first  classify  the  plots  of  White  Sands  into  different  plant  community  types  based 
on  their  species  compositions.  Then,  we  simulate  the  population  and  diversity  of 
each  type  of  plant  community  based  on  data  collected  from  the  field. 


Classification  of  Piant  Communities 

A  plant  community  is  an  interacting  unit  of  all  the  populations  of  species  within 
a  prescribed  region.  Different  communities  have  different  species  composition 
and  present  different  patterns.  A  central  goal  of  plant  ecology  and  biology  is  to 
understand  and  model  the  processes  and  mechanisms  that  cause  the  patterns  we 
see.  This  is  no  easy  task.  Any  given  habitat  may  contain  from  a  few  plant  spe¬ 
cies  to  hundreds  of  different  species.  In  a  large  area,  especially  a  tropical  or  sub¬ 
tropical  area,  species  composition  is  usually  very  complicated  and  vegetation 
patterns  are  complex.  No  single  model  can  produce  a  satisfactory  description  for 
the  population  dynamics  in  all  types  of  plant  communities.  A  model  usually 
works  well  within  a  certain  t3^e  of  community.  Therefore,  to  xmderstand  the 
vegetation  changes  in  a  given  area,  the  first  task  we  have  to  accomplish  is  to 
classify  that  community. 

In  most  sampling  designs,  vegetation  types  or  land  cover  types  are  identified 
based  on  satellite  imagery,  soil  types,  and  vegetation  information.  Samples  are 
then  proportionally  assigned  to  the  land  cover  categories  classified  from  satellite 
imagery.  For  a  vauiety  of  reasons,  classifications  based  on  satellite  imagery  do 
not  reflect  the  natural  distribution  of  plants  and  land  covers.  The  preliminary 
classification  based  on  satellite  imagery  should  be  modified  by  observation  in¬ 
formation. 

There  are  many  different  approaches  to  and  kinds  of  commxmity  classification. 
The  two  basic  methods  are  the  ecological  classification  framework  and  numerical 
classification.  The  ecological  classification  framework  combines  the  physical  and 
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biological  factors  of  a  given  region.  Usually,  dominant  species  of  a  plant  commu¬ 
nity  are  used  to  represent  the  ecological  function  of  the  commmiity.  Thus,  a  di¬ 
rect  method  for  ecological  classification  is  to  put  samples  into  groups  of  domi¬ 
nant-types  according  their  species  composition.  That  is,  samples  with  similar 
major  species  are  classified  into  the  same  group.  This  method  is  direct  and  accu¬ 
rate  when  all  samples  are  dominated  by  a  few  species.  However,  when  the  sam¬ 
ple  size  is  very  large  and  species  composition  is  very  comphcated,  the  direct 
method  may  give  a  large  number  of  dominant  types.  If  the  number  of  dominant 
types  is  limited,  it  is  very  hard,  even  for  an  experienced  expert,  to  tell  which 
group  a  sample  should  go  into.  Moreover,  this  ecological  method  is  more  or  less 
subjective.  Different  people  may  have  different  classification  schemes  for  the 
same  data  set. 

In  the  past  three  or  foiu*  decades,  there  has  been  an  increasing  tendency  to  use 
numerical  methods.  This  has  largely  been  due  to  the  greater  objectivity  of  these 
methods.  A  common  method  is  cluster  analysis.  Clustering  is  a  method  of  find¬ 
ing  groups  in  data.  But  the  direct  use  of  species  abundance  data  in  clustering 
causes  the  problem  of  putting  samples  with  similar  species  composition  into  dif¬ 
ferent  groups.  For  example,  two  samples  (I  and  II)  both  have  species  A,  B,  C, 
and  D.  Sample  I  has  the  species  abundance  100,  60,  40,  and  10  of  species  A,  B, 
C,  and  D,  respectively.  Sample  II  has  the  species  abxmdance  60,  30,  20,  and  5  of 
species  A,  B,  C,  and  D,  respectively.  Although  these  two  samples  have  similar 
species  composition  they  may  be  placed  into  different  groups  by  cluster  analysis 
because  the  absolute  difference  of  species  composition  between  these  two  sam¬ 
ples  is  large. 

Instead  of  using  the  direct  abimdance,  we  standardize  the  abundance  data  by 
quadrate  in  this  study.  Then,  we  apply  usual  clustering  methods  to  the  stan¬ 
dardized  data.  This  method  overcomes  the  disadvantages  of  using  direct  abim¬ 
dance  because  the  abundance  of  species  is  standardized  to  a  scale  of  unit  one  for 
all  samples.  The  clustering  programming  is  written  in  the  SAS  statistical  analy¬ 
sis  program.  We  use  three  different  clustering  methods:  average,  centroid,  and 
Ward’s  minimum  variance,  and  choose  the  best  clustering  scheme  after  compar¬ 
ing  the  clustering  results  obtained  by  these  three  methods. 

Standardization  of  Data  by  Quadrat 

There  are  many  different  approaches  to  standardizing  data  onto  a  scale  of  one. 
We  used  quadratic  standardization  (scaling)  in  this  study.  The  reason  for  using 
quadratic  scaling  is  that  this  method  of  scaling  puts  appropriate  weights  on 
abundant  species.  Relative  abundance,  which  is  defined  as  the  i-th  species 
abundance  divided  by  the  total  population,  treats  all  the  species  as  having  the 
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same  importsmce  in  the  classification.  This  is  not  what  we  expect  in  the  classifi¬ 
cation  of  a  plant  commmiity.  In  community  classification,  abimdant  species  are 
usually  considered  to  be  more  important  than  other  species.  On  the  contrary, 
cubic  or  higher  order  scaling  puts  too  much  weight  on  the  abimdant  species,  and 
other  species  may  have  little  or  no  impact  on  the  classification. 

To  imderstand  quadratic  scaling,  we  first  define  the  length  of  a  vector.  The 
length  of  a  vector  is  defined  as  the  square  root  of  the  sum  of  the  squares  of  all  its 
elements.  That  is, 

(1) 

where  is  the  i-th  element  of  vector  Q. 

Let  Q  be  the  abimdance  vector  of  a  community.  is  the  abxmdaince  of  i-th  spe¬ 
cies.  Then,  the  standardized  abxmdance  of  the  i-th  species,  which  is  also  called 
the  importance  index  of  the  i-th  species,  is  defined  by: 

(2)  Qi  =  jrQf 

Two  properties  of  the  importance  index  are:  (1)  Q'  is  a  unit  vector;  that  is, 
ranges  from  0  to  1,  and  =  1;  (2)  is  determined  by  the  relative  abimdance 
rather  than  absolute  abundance. 

Clustering  Method 

The  SAS  clustering  procedure,  CLUSTER,  is  used  to  find  groups  of  observations 
with  coordinate  data  (species  importance  index).  To  obtain  a  better  clustering 
result,  we  use  three  clustering  methods:  average  linkage,  centroid  method,  and 
Ward’s  minimum-variance  method. 

The  following  notation  is  used,  with  lowercase  S3rmbols  generally  pertaining  to 
observations  and  uppercase  s3rmbols  to  clusters: 
n  =  number  of  samples  (observations); 

V  =  number  of  variables; 

G  =  number  of  cliisters  at  any  given  level  of  the  hierarchy; 
x^  =  i-th  observation; 

=  K-th  cluster; 

=  nimiber  of  observations  in  C^; 

Xk  =  mean  vector  for  C^; 
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X  =  Euclidean  length  of  the  vector  x; 


=  E/6C.  p/ -  xa  I  ; 

=Wj^-W^-W,ifC^=Q^C,; 

D(x,  y)  =  distance  between  vectors  x  and  y; 

=  distance  between  clusters  and 

The  distance  between  two  clusters  can  be  defined  either  directly  or  by  combina¬ 
torial.  That  is,  by  an  equation  for  updating  a  distance  matrix  when  two  clusters 
are  joined.  In  all  combinatorial  formulas  below,  it  is  assumed  that  clusters 
and  Cl  are  merged  to  form  C^,  and  the  formula  gives  the  distance  between  the 
new  cluster  C„  and  any  other  cluster  Cj. 

Average  Linkage 

In  the  average  linkage,  the  distance  between  two  clusters  is  defined  by: 

(3)  Dfci  —  (NkNl). 

if  d(x,y)  =  \\x  -  yf  then: 

(4)  Dkl  =  Xk-Xi  +Wk  I  'Nk'^Wi  I  Ni. 

The  combination  formula  is: 

(5) 

In  the  average  linkage,  the  distance  between  two  clusters  is  the  average  distance 
between  pairs  of  observations,  one  in  each  cluster.  Average  hnkage  tends  to  join 
clusters  with  small  variances  but  is  slightly  biased  toward  producing  clusters 
with  the  same  variance. 

Centroid  Method 

In  the  centroid  method,  the  distance  between  two  clusters  is  defined  by: 

Dkl  ~  “  Xi  I  . 


(6) 
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If  d(x,y)  =  ||x  -  vir  combination  formula  is: 

(7)  Djm  =  {J^kDjK  +  )I'Nm-  ^k'NlDkl  /  • 

In  the  centroid  method  the  distance  between  two  clusters  is  defined  as  the 
Euclidean  distance  between  their  centroid  or  means.  The  centroid  method  is 
more  robust  to  outliers  than  most  other  hierarchical  methods. 

Ward’s  Minimum-variance  Method 

In  Ward’s  minimum-variance  method,  the  distance  between  two  clusters  is  the 
analysis  of  variance  (ANOVA)  sum  of  squares  between  two  clusters  added  up 
over  all  the  variables.  At  each  generation,  the  within-cluster  sum  of  squares  is 
minimized  over  all  partitions  obtainable  by  merging  two  clusters  from  the  previ¬ 
ous  generation.  d(x,  y),  and  the  combinatorial  formulas  are,  respectively,  de¬ 
fined  as: 

(8)  d(x,y)=  I  Ix-yl  IV2 

(9)  D^=  I  IXk-XlI  iVd/N^+l/N,) 

(10)  =  ((N,  -H  -H  (N,  -h  -  N,D^)  /  (N,  -h  N^) 

Classification  Results 

The  clustering  results  were  obtained  by  the  three  clustering  methods:  average 
linkage,  centroid  method,  and  Ward’s  minimum-variance  method.  Tables  1  and  2 
summarize  the  results  of  cluster  analysis  for  years  1989  and  1992,  respectively. 
From  Tables  1  and  2  we  can  see  Ward’s  minimiun  variance  method  gives  the  best 
match  of  plots  among  all  the  three  methods  and  a  stable  classification  scheme 
between  years  1989  and  1992.  Thus  we  choose  Ward’s  clustering  scheme  as  the 
frame  scheme  and  adjusted  it  by  the  clustering  results  from  average  linkage  and 
centroid  method.  We  suggest  18  vegetation  types  for  White  Sands  (17  clusters 
plus  1  type  of  bare  soil).  The  vegetation  types  and  their  major  species  are  listed 
in  Table  3. 
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Cluster  Average  Linkage 


1  31,  33,  35,  36,  147,  150,  175 


2  1,  2,  10,  11,  12,  13,  14,  24, 

29,  30,  37,  41,  42,  54,  59,  73, 
80,  82,  83.  84,  85,  86,  88.  89, 
90,  92,  93,  95,  96,  110,  113, 
117,  118,  119,  120,  121,  122, 
123,  125,  132,  133,  134,  151, 
156,  157,  158,  159,  161,  166, 
180,  181,  189,  190,  191,  196 


3  19,  20,  26,  53,  57,  58,  68,  69, 

91,  109,  116,  129,  131,  162, 

168,  176,  200 _ 

3.  4,  23,  25,  27,  38,  39.  56, 

60,  61,  62,  114,  115,  127,  152, 
155,  172,  195 


5  46,  49,  100,  101,  102,  112, 

177,  197 


6  6,  7,  8,  21,  22,  51,  52,  64, 

128,  130,  138,  139,  144,  160, 

164,  165,  167 _ 

70,  111,  142,  145,  171,  173, 

199 _ 

8  15,  16,  94,  97,  124,  135,  178 


Centroid  Hierarchical 


131,  33,  35,  36.  147,  150,  175 


1,  2,  10.  11,  12,  13,  14,  24, 
29,  30,  41,  42,  54,  59,  73,  82, 
83,  84,  85.  86,  88,  89,  90,  92, 
93.  110,  113,  117,  118,  119, 
120,  121,  125,  132,  133,  151, 
156,  157,  161,  166,  180,  181, 
189,  191 


Ward’s  Minimum  Variance 


|31,  33,  35,  36.  147,  150,  175 


1,  10,  12,  14,  24,  29,  30,  41, 
42,  54,  59,  82.  84.  85,  86,  88. 
89.  90,  92,  93,  113,  117,  118, 
119,  120,  121,  125,  132,  133, 
156,  157,  161,  166,  189,  191 


153,  57,  68,  69,  91,  109,  129  |53,  57,  68,  69,  91,  109,  129 


5,  9,  65,  146,  153,  154,  193, 
194 


3,  4.  23,  25,  27,  38,  39,  56, 

60,  62,  114,  115,  127,  152, 

155,  172,  195  _ 


46,  49.  100,  101,  102,  112, 

177,  197  _ 


7,  21.  22,  64,  138,  139,  144, 
160,  164,  165,  167 

70,  111,  142,  145,  171,  173, 

199 _ 

6,  8,  15,  16,  19,  20,  26,  37, 

47,  50,  51,  52,  58,  61,  66,  67, 

71,  72,  75,  77,  78,  79,  80,  87, 

94,  95,  96,  97,  99,  106,  116, 
122,  123,  124,  128,  130,  131, 

134,  135,  137,  158,  159,  162, 

163,  168,  176,  178,  182,  184, 

190,  196,  198,  200 _ 

5,  9,  65,  146,  153,  154 


10  40,  48.  79,  98,  103,  104,  108, 

183 


11  105,  107  105,  107 


12  28,  43,  126 


193,  194 


47,  50,  66,  67,  71,  72,  75,  77,  28,  43,  126 
78,  87,  99,  106,  137,  163,  182, 

184,  198 


17,  18,  81 


17,  18,  81 


74 _ 

63 


'6 _ 

36 


140 


169 


3,  4,  23,  25,  27,  38,  39,  56, 

60,  62,  114,  115,  127,  152, 

155,  172,  195  _ 


46,  49,  100,  101,  102,  112, 

177,  197  _ 


21,  22,  64,  138,  139,  160,  164, 
167 


70,  111,  142,  145,  171,  173, 

199 _ 

19,  20,  26,  58,  116,  131,  162, 
176,  200 


2,  11,  13,  37,  73.  83,  110,  122, 
123,  134,  151,  158,  180,  181, 
190 


5,  £ 

65.  146,  153, 

154, 

193, 

194 

6,  7 

8. 

51, 

52, 

61, 

128, 

130, 

144 

159,  165, 

196 

50, 

63. 

66. 

67. 

71, 

72,  74,  75, 

76, 

77, 

78, 

136 

137,  163,  182, 

184 

105 

107 

28, 

43, 

126 

15, 

16, 

94, 

97. 

124, 

135, 

178 

17, 

18, 

81 

47, 

87, 

99, 

106, 

198 

140 

169 
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Table  2.  Comparison  of  cluster  analysis  results  for  1992. 


Cluster 


Average  Linkage 

Centroid  Hierarchical 

31,  33.,  36,  147,  148,  175 

31,  33,  36,  147,  148,  175 

1,  10.  11,  12,  13,  14,  24,  29, 

1,  2,  10,  11,  12,  13,  14,  24, 

30,  41,  42,  54.  59.  73,  82,  88, 

29,  30,  41.  42,  54,  59,  73,  82, 

89,  90,  92,  93,  113,  117,  118, 

88,  89,  90,  92,  93,  113,  117, 

119,  120,  121,  123,  125,  132, 

118,  119,  120,  121,  123,  125, 

133,  156,  157,  158,  161,  166, 

132,  133,  156,  157,  158,  161, 

180,  181,  189,  191 

166,  180,  181,  189,  191 

3,  4,  19,  20,  23,  25,  26,  27. 

3,  4,  23,  27,  38,  39,  60,  62, 

38,  39,  56,  57,  58.  60.  62,  68, 
69,  91,  109,  114,  115,  116, 

127,  129,  130,  131,  152,  162, 
168,  172,  195,  200 

127,  152,  172,  195  . 

70,  77,  111,  142,  145,  151, 

70,  77,  111,  142,  145,  151, 

171,  173,  199 

171,  173,  199 

46,  49.  100,  101,  102,  110, 

46,  49,  100,  101,  102,  110, 

112,  177,  197 

112,  177,  197 

6,  7,  17,  21,  22,  35,  47, 

6,  7,  15,  16,  17.  19,  20,  21, 

50.  51,  61,  64,  66,  67, 

22,  25,  26,  37,  40,  48.  50,  51, 

71,  72,  74,  75,  78,  99, 

56,  57.  58,  64,  66,  67,  68,  69, 

71,  72,  74,  75,  78,  79.  80,  91, 

103,  106,  128,  137,  138, 

94,  95,  96,  97,  98,  104,  105, 

139,  140,  141,  144,  159, 

108,  109,  114,  115,  116,  122, 

160,  163,  164,  165,  167, 

124,  128,  129,  130,  131,  134, 

182,  183,  184,  196,  198 

135,  137,  138,  139,  140,  141, 
144,  159,  160,  162,  163,  164, 
165,  167,  168,  176,  178,  182, 
183,  184,  190,  196,  200 

2,  15,  16,  37,  40,  48,  79,  80, 
94,  95,  96,  97,  98.  104,  105, 
108,  122,  124,  134,  135,  176, 
178,  190 

35,  47,  61,  99,  103,  106,  198 

5,  65.  146,  153 

5,  65,  146,  153 

28,  43,  126 

28.  43,  126 

193,  194 

193,  194 

63 

63 

107 

107 

76 

76 

18 

18 

81 

81 

136 

136 

8 

8 

154 

154 

169 

169 

52 

52 

Ward’s  Minimum  Variance 


31,  33,  36,  147,  148,  175 


1,  10,  11,  12,  13,  14,  24,  29, 
30,  41.  42,  54,  59,  73,  82,  88, 
89,  90,  92,  93,  113,  117,  118, 

119,  120,  121,  123,  125,  132, 

133,  156,  157,  158,  161,  166, 

180,  181,  189,  191 


3,  4,  23,  27,  38,  39,  60,  62, 
127,  152,  172,  195 


70,  77,  111,  142,  145,  151, 
171,  173,  199 


46,  49,  100,  101,  102,  110, 
112,  177,  197  _ 


6,  7,  21,  22.  51,  64,  128, 
138,  139,  159,  160,  164, 
165,  167,  184,  196 


15,  16,  37,  94,  95,  96.  97,  122, 
124,  134,  135,  176,  178,  190 


26,  57,  58,  68,  69,  91,  109, 
129,  162,  168,  200 


19,  20,  25,  56,  114,  115,  116, 
130,  131 


2,  40,  48.  79,  80,  98,  104,  105, 
108 


5,  65,  146,  153 


28,  43.  126 


35,  47,  61,  99,  103,  106,  183, 
198 


193,  194 


50,  63,  66.  67,  72,  74,  75,  76, 
78.  107,  136,  137,  140,  141, 
163,  182 


17,  18.  71,  81,  144 
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Table  3.  Vegetation  types  identified  at  White  Sands  Missile  Range. 


IDB 

No.  of  piots 

Piot  ID 

Major  species 

1 

3 

36,  147,  175 

ALOC2,  PSAR,  ATCA2 

2 

35 

1,  10,  12,  14,  24,  29,  30,  41,  42,  54,  59, 

82,  84,  85,  86,  88,  89,  90,  92,  93,  113,  117, 
118,  119,  121,  125,  132,  133,  156,  157,  161, 
166,  189,  191 

ATCA2 

3 

7 

53,  57,  68,  69,  91,  109,  129 

FLCE 

H 

17 

3,  4,  23,  25,  27,  38,  39,  56,  60,  62,  114, 

115,  127,  152,  155,  172,  195 

LATR2 

8 

46,  49,  100,  101,  102,  112,  177,  197 

ARFI2,  EPTO 

6 

8 

21,  22,  64,  138,  139,  160,  164,  167 

PAIN2,  LATR2,  VIST.  DAFO 

7 

7 

70,  111,  142,  145,  171,  173,  199 

DAFO,  PAIN2,  LATR2 

8 

9 

19,  20,  26,  58,  116,  131,  162,  176,  200 

FLCE,  LATR2 

9 

14 

2,  11,  13,  37,  73,  83,  122,  123,  134,  151, 

158,  180,  181,  190 

ATCA2,  LYBE 

10 

6 

15,  94,  97,  124,  135,  178 

LYBE,  ATCA2 

11 

8 

5,  9,  65,  146,  153,  154,  193,  194 

VIST,  DAFO,  ACC02,  PAIN2 

12 

12 

6,  7,  8,  51,  52,  61,  128,  130,  144,  159,  165, 
196 

PAIN2.  LATR2,  DAFO 

13 

12 

40,  48,  79,  80,  95,  96,  98,  103,  104,  108, 

168,  183 

EPTO,  ATCA2 

B 

18 

50,  63,  66,  67,  71,  74,  75,  76,  77,  78,  105, 
107,  136,  137,  140,  163,  182,  184 

DAWH2,  PAIN2,  YUBA 

mm 

3 

17,  18,  81 

OPVI,  ACGR 

4 

47,  87,  106,  198 

YUEL,  POIN3 

17 

1 

169 

DAFO.  PAIN2,  FLCE 

18 

27 

16,  28,  31,  33,  35,  43,  72,  99,  110, 120, 126, 150 
32,  34,  44,  45,  55, 143,  148,  149,  174, 179,  185, 
186, 187, 188, 192 

Bare  land 
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Simulation  Results 


The  STELLA  model  as  described  in  Chapter  5  was  used  to  model  the  population 
and  diversity  of  the  17  plant  community  types.  The  data  used  for  the  simulation 
were  from  four  years:  1989,  1990,  1991,  and  1992.  Figures  4  through  20  illus¬ 
trate  results  of  the  simulation  for  the  different  plant  communities. 


1:  TOTAL  POPULATION 


2:  TOTAL  SPECIES 


3:  SHANNON  INDEX 


M: 
:  2: 
P3: 


i  2: 
.3: 


200.00- 

20.00 

3.00 


1(X).00 

10.00- 

1.50 

■■ 


0.00 

0,00 

0.00- 


,  '■  *  f''r\  .  \-,<i 

!'  r  ^  '' ■  •. 

!4.r- . ■  • " r  s' '*»=:■  1..V  ’ ■ 

•  ■  .  ■  .'Ml';..;-  .  ;■  •  ■ 

d..--  ... 

2  :  . 

0.00 

0  Graph:  Page  1 


4.00 


8.00 

Years 


12,00  16.00 
11:49  AM  11/29/95 


Figure  5.  Simulation  results  for  plant  community  type  2. 
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1:  TOTAL  POPULATION  2;  TOTAL  SPECIES  3;  SHANNON  INDEX 


Figure  6.  Simulation  results  for  plant  community  type  3. 
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Figure  8.  Simulation  results  for  plant  community  type  5. 
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Figure  9.  Simulation  results  for  plant  community  type  6. 
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1:  TOTAL  POPULATION 
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0  Graph:  Page  1 
Figure  10.  Simulation  results  for  plant  community  type  7 
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Figure  11.  Simulation  results  for  plant  community  type  8. 
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1:  TOTAL  POPULATION  2:  TOTAL  SPECIES  3:  SHANNON  INDEX 


Figure  12.  Simulation  results  for  plant  community  type  9. 
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Figure  16.  Simulation  results  for  plant  community  type  13. 


Figure  17.  Simulation  results  for  plant  community  type  14. 
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Figure  18.  Simulation  results  for  plant  community  type  15. 
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Figure  20.  Simulation  results  for  plant  community  type  17. 
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7  Summary 

This  report  is  the  cuhnination  of  a  project  that  was  designed  to  develop  and  test 
a  new  methodology  to  model  changes  in  plant  diversity.  Using  standard  data 
from  the  U.S.  Army’s  LCTA  program  at  White  Sands  Missile  Range,  New  Mexico, 
stochastic  models  of  plant  diversity  were  created  to  simulate  the  dynamics  of  the 
population  growth  in  a  nximber  of  plant  communities  on  the  installation. 

This  project  resulted  in  the  development  of  a  new  model  in  which  we  demon¬ 
strated  different  methods  to  solve  a  variety  of  stochastic  differential  equations. 
We  illustrated  how  to  incorporate  both  demographic  noise  and  environmental 
noise  into  a  single  model  containing  the  joint  effects  of  demographic  and  envi¬ 
ronmental  stochasticity.  This  new  model  contains  two  main  components:  deter¬ 
ministic  process  and  stochastic  process.  We  first  introduced  a  dynamics  model  of 
population  growth.  This  model  derives  birth  and  death  rates  as  they  relate  to 
population  growth  from  the  relationships  among  plemts.  It  also  derives  these 
rates  from  the  relationship  between  plants  and  the  environment.  This  model 
serves  as  the  deterministic  part  of  the  stochastic  model  and  controls  the  biologi¬ 
cal  d3Taamics  of  the  plant  communities.  The  second  component  of  the  model 
simulates  the  biological  and  ecological  fluctuations.  In  the  new  model,  the  sto¬ 
chastic  process  is  simulated  with  the  birth  and  death  process.  This  process  best 
describes  the  demographic  stochasticity  because  it  depicts  popxilation  dynamics 
with  biologically  accurate,  interpretable  birth  and  death  rates  and  is  apphcable 
to  populations  of  var3dng  sizes. 

Diiring  the  testing  phase  of  the  project,  we  needed  to  have  clearly  defined  plant 
commiuiities  in  order  to  test  the  model  with  White  Sands  LCTA  data.  We  used 
Ward’s  minimum  variance  method  of  cluster  analysis  because  it  yielded  the  best 
match  of  plots  out  of  the  three  methods  of  cluster  analysis  used  to  characterize 
individual  plant  communities  on  the  White  Sands  installation. 

This  diversity  model  complements  other  models  that  use  LCTA  data  to  deter¬ 
mine  plant  population  levels.  Plant  population  models  are  useful  tools  in  that 
they  enable  natural  resource  managers  to  transcend  current  identification 
strategies  and  can  help  determine  future  training  levels  that  will  allow  the 
maximum  level  of  training  to  occur  in  areas  with  minimal  impact  on  overall  spe¬ 
cies  diversity.  Although  this  particular  model  is  only  applicable  at  the  White 
Sands  installation,  this  new  model  can  handle  a  variable  niunber  of  species  and 
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species  abundance.  This  capability  enables  it  to  be  used  to  calculate  any  kind  of 
diversity  indices.  Because  both  the  dynamic  changes  and  stochastic  fluctuations 
are  included  in  the  stochastic  dynamics  model,  re-parameterizing  the  model 
gives  it  the  potential  for  extensive  use  in  natural  resource  management  and  en¬ 
vironmental  monitoring.  Future  applications  of  these  stochastic  dynamics  mod¬ 
els  include:  providing  standard  diversity  measimes;  monitoring  the  development 
of  plant  communities  in  terms  of  species  diversity  and  structure  diversity;  test¬ 
ing  the  significance  of  the  influence  of  human  activities  on  plant  commumties; 
and  estimating  rehabilitation  time  for  a  distxirbed  plant  conmnmity,  thereby 
helping  the  Army  integrate  training  and  natural  resource  management. 
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